2.2 Gamma Function

Gamma function: generalization of factorial function to the real or complex numbers. Factorial numbers are only defined with 0 or positive integers. Here we only handle real number Gamma function. For real work, use gamma R internal function. Many distribution and probability functions depend on gamma or log gamma function.

\[\Gamma (x) = \int_0^{\infty} t^{x-1}\mathrm{e}^{-t}\mathrm{d}t\]

while x is from 0 to infinity.

See Rudin’s book 3e p192 for more detail.

User defined GAMMA function using Lanczos approximation

GAMMA

## function (z) 
## {
##     if (is.nan(z)) 
##         return(NaN)
##     if (z == -Inf) 
##         return(NaN)
##     if (z > 171.61) 
##         return(+Inf)
##     if (z <= 0 & z == floor(z)) 
##         return(NaN)
##     if (z < 0.5) 
##         return(pi/(sin(pi * z) * GAMMA(1 - z)))
##     if (z == floor(z) & z < 172) 
##         return(tableFactorial(z - 1))
##     sqrt2pi = 2.506628274631
##     p = c(676.520368121885, -1259.1392167224, 771.323428777653, 
##         -176.615029162141, 12.5073432786869, -0.13857109526572, 
##         9.98436957801957e-06, 1.50563273514931e-07)
##     z = z - 1
##     x = 0.99999999999981
##     for (i in 1:8) x = x + p[i]/(z + i)
##     t = z + 7.5
##     if (z > 141.2) {
##         return(sqrt2pi * exp((z + 0.5) * log(t) - t + log(x)))
##     }
##     else {
##         return(sqrt2pi * t^(z + 0.5) * exp(-t) * x)
##     }
## }
## <bytecode: 0x000001e2a88f6540>
## <environment: namespace:math>

GAMMA(0)

## [1] NaN

GAMMA(1)   # 0! = 1

## [1] 1

GAMMA(2)   # 1! = 1

## [1] 1

GAMMA(171) # 170! = 7.257416e+306

## [1] 7.257e+306

format(GAMMA(0.5), digits=22) # = sqrt(pi) = 1.772453850905516

## [1] "1.772453850905516326009"

format(gamma(0.5), digits=22)

## [1] "1.772453850905516103964"

format(sqrt(pi), digits=22)

## [1] "1.772453850905515881919"

format(GAMMA(1.5), digits=22)

## [1] "0.8862269254527587181158"

format(gamma(1.5), digits=22)

## [1] "0.886226925452758051982"

For very large values (the result value exceeding .Machine$double.xmax), we need to use log gamma function.

LGAMMA

## function (z) 
## {
##     if (is.nan(z)) 
##         return(NaN)
##     if (z == -Inf) 
##         return(NaN)
##     if (z > 2.53273727608008e+305) 
##         return(+Inf)
##     if (z <= 0 & z == floor(z)) 
##         return(NaN)
##     if (z < 0.5) 
##         return(log(pi/sin(pi * z)) + LGAMMA(1 - z))
##     if (abs(z) <= 171) 
##         return(log(abs(GAMMA(z))))
##     else if (abs(z) > 1e+17) 
##         return(z * (log(z) - 1))
##     else {
##         lnsqrt2pi = 0.918938533204673
##         p = c(676.520368121885, -1259.1392167224, 771.323428777653, 
##             -176.615029162141, 12.5073432786869, -0.13857109526572, 
##             9.98436957801957e-06, 1.50563273514931e-07)
##         z = z - 1
##         x = 0.99999999999981
##         for (i in 1:8) x = x + p[i]/(z + i)
##         t = z + 7.5
##         return(lnsqrt2pi + (z + 0.5) * log(t) - t + log(x))
##     }
## }
## <bytecode: 0x000001e2a3ffc9f8>
## <environment: namespace:math>

LGAMMA(0)

## [1] NaN

LGAMMA(1)   # log(0!) = 0

## [1] 0

LGAMMA(2)   # log(1!) = 0

## [1] 0

LGAMMA(171) # log(170!) = 706.5731

## [1] 706.6

LGAMMA(-171) # NaN for all negative intergers

## [1] NaN

format(LGAMMA(0.5), digits=22) # = log(sqrt(pi)) = 0.5723649429247

## [1] "0.5723649429247003039833"

format(lgamma(0.5), digits=22)

## [1] "0.5723649429247000819387"

format(log(sqrt(pi)), digits=22)

## [1] "0.5723649429246999709164"

format(LGAMMA(1.5), digits=22)

## [1] "-0.1207822376352444271319"

format(lgamma(1.5), digits=22)

## [1] "-0.1207822376352451765325"

When you call factorial function in R, it calls gamma function.

Comparison of LGAMMA function and lgamma R internal funciton

Start = 3      # skip zero values
End = 1000000
RelTol = 1e-15 # relative tolerance
for (i in Start:End) {
  x1 = LGAMMA(i)
  x2 = lgamma(i)
  RelDiff = abs((x1 - x2)/x2)
  if (RelDiff > RelTol) print(paste(i, format(c(x1, x2), digits=22)))
}

AbsTol = 1e-8 # absolute tolerance
for (i in Start:End) {
  x1 <- LGAMMA(i)
  x2 <- lgamma(i)
  if (abs(x1 - x2) > AbsTol) print(paste(i, format(c(x1, x2), digits=22)))
}

With above test, there should be no printout to be successful. If you make the tolerance more strict, you can see the differences.