Special Functions of Mathematics

Description

Special mathematical functions related to the beta and gammafunctions.

Usage

beta(a, b)lbeta(a, b)gamma(x)lgamma(x)psigamma(x, deriv = 0)digamma(x)trigamma(x)choose(n, k)lchoose(n, k)factorial(x)lfactorial(x)

Arguments

Details

The functionsbeta andlbeta return the beta functionand the natural logarithm of the beta function,

B(a,b) = Γ(a)Γ(b)/Γ(a+b).

The formal definition is

integral_0^1 t^(a-1) (1-t)^(b-1) dt

(Abramowitz and Stegun section 6.2.1, page 258). Note that it is onlydefined inR for non-negativea andb, and is infiniteif either is zero.

The functionsgamma andlgamma return the gamma functionΓ(x) and the natural logarithm ofthe absolute value of thegamma function. The gamma function is defined by(Abramowitz and Stegun section 6.1.1, page 255)

Γ(x) = integral_0^Inf t^(x-1) exp(-t) dt

for all realx except zero and negative integers (whenNaN is returned). There will be a warning on possible loss ofprecision for values which are too close (within about1e-8) to a negative integer less than-10.

factorial(x) (x! for non-negative integerx)is defined to begamma(x+1) andlfactorial to belgamma(x+1).

The functionsdigamma andtrigamma return the first and secondderivatives of the logarithm of the gamma function.psigamma(x, deriv) (deriv >= 0) computes thederiv-th derivative ofψ(x).

digamma(x) = ψ(x) = d/dx{ln Γ(x)} = Γ'(x) / Γ(x)

ψ and its derivatives, thepsigamma() functions, areoften called the ‘polygamma’ functions, e.g. inAbramowitz and Stegun (section 6.4.1, page 260); and higherderivatives (deriv = 2:4) have occasionally been called‘tetragamma’, ‘pentagamma’, and ‘hexagamma’.

The functionschoose andlchoose return binomialcoefficients and the logarithms of their absolute values. Note thatchoose(n, k) is defined for all real numbersn and integerk. Fork ≥ 1 it is defined asn(n-1)…(n-k+1) / k!,as1 fork = 0 and as0 for negativek.Non-integer values ofk are rounded to an integer, with a warning.
choose(*, k) uses direct arithmetic (instead of[l]gamma calls) for smallk, for speed and accuracyreasons. Note the functioncombn (packageutils) for enumeration of all possible combinations.

Thegamma,lgamma,digamma andtrigammafunctions are internal generic primitive functions: methods can bedefined for them individually or via theMath group generic.

Source

gamma,lgamma,beta andlbeta are based onC translations of Fortran subroutines by W. Fullerton of Los AlamosScientific Laboratory (now available as part of SLATEC).

digamma,trigamma andpsigamma are based on

Amos, D. E. (1983). A portable Fortran subroutine forderivatives of the psi function, Algorithm 610,ACM Transactions on Mathematical Software9(4), 494–502.

References

Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988)The New S Language.Wadsworth & Brooks/Cole. (Forgamma andlgamma.)

Abramowitz, M. and Stegun, I. A. (1972)Handbook of Mathematical Functions. New York: Dover.https://en.wikipedia.org/wiki/Abramowitz_and_Stegun provideslinks to the full text which is in public domain.
Chapter 6: Gamma and Related Functions.

See Also

Arithmetic for simple,sqrt formiscellaneous mathematical functions andBessel for thereal Bessel functions.

For the incomplete gamma function seepgamma.

Examples

require(graphics)choose(5, 2)for (n in 0:10) print(choose(n, k = 0:n))factorial(100)lfactorial(10000)## gamma has 1st order poles at 0, -1, -2, ...## this will generate loss of precision warnings, so turn offop <- options("warn")options(warn = -1)x <- sort(c(seq(-3, 4, length.out = 201), outer(0:-3, (-1:1)*1e-6, "+")))plot(x, gamma(x), ylim = c(-20,20), col = "red", type = "l", lwd = 2,     main = expression(Gamma(x)))abline(h = 0, v = -3:0, lty = 3, col = "midnightblue")options(op)x <- seq(0.1, 4, length.out = 201); dx <- diff(x)[1]par(mfrow = c(2, 3))for (ch in c("", "l","di","tri","tetra","penta")) {  is.deriv <- nchar(ch) >= 2  nm <- paste0(ch, "gamma")  if (is.deriv) {    dy <- diff(y) / dx # finite difference    der <- which(ch == c("di","tri","tetra","penta")) - 1    nm2 <- paste0("psigamma(*, deriv = ", der,")")    nm  <- if(der >= 2) nm2 else paste(nm, nm2, sep = " ==\n")    y <- psigamma(x, deriv = der)  } else {    y <- get(nm)(x)  }  plot(x, y, type = "l", main = nm, col = "red")  abline(h = 0, col = "lightgray")  if (is.deriv) lines(x[-1], dy, col = "blue", lty = 2)}par(mfrow = c(1, 1))## "Extended" Pascal triangle:fN <- function(n) formatC(n, width=2)for (n in -4:10) {    cat(fN(n),":", fN(choose(n, k = -2:max(3, n+2))))    cat("\n")}## R code version of choose()  [simplistic; warning for k < 0]:mychoose <- function(r, k)    ifelse(k <= 0, (k == 0),           sapply(k, function(k) prod(r:(r-k+1))) / factorial(k))k <- -1:6cbind(k = k, choose(1/2, k), mychoose(1/2, k))## Binomial theorem for n = 1/2 ;## sqrt(1+x) = (1+x)^(1/2) = sum_{k=0}^Inf  choose(1/2, k) * x^k :k <- 0:10 # 10 is sufficient for ~ 9 digit precision:sqrt(1.25)sum(choose(1/2, k)* .25^k)