# Gamma and Beta function

### From HaskellWiki

(Difference between revisions)

m |
|||

(5 intermediate revisions by 2 users not shown) | |||

Line 2: | Line 2: | ||

<haskell> |
<haskell> |
||

cof :: [Double] |
cof :: [Double] |
||

− | cof = [76.18009172947146,-86.50532032941677,24.01409824083091,-1.231739572450155,0.01208650973866179,-0.00005395239384953] |
+ | cof = [76.18009172947146,-86.50532032941677,24.01409824083091,-1.231739572450155,0.001208650973866179,-0.000005395239384953] |

ser :: Double |
ser :: Double |
||

Line 9: | Line 9: | ||

gammaln :: Double -> Double |
gammaln :: Double -> Double |
||

gammaln xx = let tmp' = (xx+5.5) - (xx+0.5)*log(xx+5.5) |
gammaln xx = let tmp' = (xx+5.5) - (xx+0.5)*log(xx+5.5) |
||

− | ser' = foldl (+) ser $ map (\(y,c) -> c/(xx+y)) $ zip [1..] cof |
+ | ser' = ser + sum $ zipWith (/) cof [xx+1..] |

− | in -tmp' + log(2.5066282746310005 * ser' / xx) where |
+ | in -tmp' + log(2.5066282746310005 * ser' / xx) |

</haskell> |
</haskell> |
||

− | the beta function relates to the gamma function by <math>B(z,w)= \frac{Gamma(z)*Gamma(w)}{Gamma(z+w)}</math>, so we can compute the Beta function using gammaln like this: |
+ | the beta function relates to the gamma function by <math>\Beta(z,w)= \frac{\Gamma(z)\cdot\Gamma(w)}{\Gamma(z+w)}</math>, so we can compute the Beta function using gammaln like this: |

<haskell> |
<haskell> |
||

− | beta z w = exp ((gammaln z) + (gammaln w) - (gammaln (z+w)) |
+ | beta z w = exp (gammaln z + gammaln w - gammaln (z+w)) |

</haskell> |
</haskell> |
||

## Latest revision as of 10:18, 13 December 2009

The Gamma and Beta function as described in 'Numerical Recipes in C++', the approximation is taken from [Lanczos, C. 1964 SIAM Journal on Numerical Analysis, ser. B, vol. 1, pp. 86-96]

cof :: [Double] cof = [76.18009172947146,-86.50532032941677,24.01409824083091,-1.231739572450155,0.001208650973866179,-0.000005395239384953] ser :: Double ser = 1.000000000190015 gammaln :: Double -> Double gammaln xx = let tmp' = (xx+5.5) - (xx+0.5)*log(xx+5.5) ser' = ser + sum $ zipWith (/) cof [xx+1..] in -tmp' + log(2.5066282746310005 * ser' / xx)

the beta function relates to the gamma function by , so we can compute the Beta function using gammaln like this:

beta z w = exp (gammaln z + gammaln w - gammaln (z+w))