Special {base} | R Documentation

## Special Functions of Mathematics

### Description

Special mathematical functions related to the beta and gamma functions.

### 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

`a, b, x, n` | numeric vectors.
`k, deriv` | integer vectors.

### Details

The functions `beta` and `lbeta` return the beta function and the natural logarithm of the beta function,

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

The formal definition is

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

(Abramowitz and Stegun (6.2.1), page 258).

The functions `gamma` and `lgamma` return the gamma function Γ(x) and the natural logarithm of the absolute value of the gamma function. The gamma function is defined by (Abramowitz and Stegun (6.1.1), page 255)

integral_0^Inf t^(a-1) exp(-t) dt

`factorial(x)` is x! and identical to `gamma(x+1)` and `lfactorial` is `lgamma(x+1)`.

The functions `digamma` and `trigamma` return the first and second derivatives of the logarithm of the gamma function. `psigamma(x, deriv)` (`deriv >= 0`) is more generally computing the `deriv`-th derivative of psi(x).

`digamma(x)` = psi(x) = d/dx {ln Gamma(x)} = Gamma'(x) / Gamma(x)

The functions `choose` and `lchoose` return binomial coefficients and their logarithms. Note that `choose(n,k)` is defined for all real numbers n and integer k. For k >= 1 as n(n-1)...(n-k+1) / k!, as 1 for k = 0 and as 0 for negative k.
`choose(*,k)` uses direct arithmetic (instead of `[l]gamma` calls) for small `k`, for speed and accuracy reasons. Note the function `combn` (package utils) for enumeration of all possible combinations.

The `gamma`, `lgamma`, `digamma` and `trigamma` functions are generic: methods can be defined for them individually or via the `Math` group generic.

### References

Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988) The New S Language. Wadsworth & Brooks/Cole. (for `gamma` and `lgamma`.)

Abramowitz, M. and Stegun, I. A. (1972) Handbook of Mathematical Functions. New York: Dover. Chapter 6: Gamma and Related Functions.

### See Also

`Arithmetic` for simple, `sqrt` for miscellaneous mathematical functions and `Bessel` for the real Bessel functions. Note that `gammaCody(x)` is considerably faster than `gamma(x)` but slightly less accurate and (potentially) less reliable.

For the incomplete gamma function see `pgamma`.

### 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 off
    op <- options("warn")
    options(warn = -1)
    x <- sort(c(seq(-3,4, length=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(.1, 4, length = 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 <- paste(ch, "gamma", sep = "")
      if (is.deriv) {
        dy <- diff(y) / dx # finite difference
        der <- which(ch == c("di","tri","tetra","penta")) - 1
        nm2 <- paste("psigamma(*, deriv = ", der,")",sep='')
        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))), "\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:6
    cbind(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)

[Package base version 2.6.0 Index]