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
non-negative numeric vectors.
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) = Γ(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 only
defined in R for non-negative a and b, and is infinite
if either is zero.
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 section 6.1.1, page 255)
Γ(x) = integral_0^Inf t^(x-1) exp(-t) dt
for all real x except zero and negative integers (when
NaN is returned). There will be a warning on possible loss of
precision for values which are too close (within about
1e-8)) to a negative integer less than -10.
factorial(x) (x! for non-negative integer x)
is defined to be gamma(x+1) and lfactorial to be
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) computes the
deriv-th derivative of ψ(x).
digamma(x) = ψ(x) = d/dx{ln Γ(x)} = Γ'(x) / Γ(x)
ψ and its derivatives, the psigamma() functions, are
often called the ‘polygamma’ functions, e.g. in
Abramowitz and Stegun (section 6.4.1, page 260); and higher
derivatives (deriv = 2:4) have occasionally been called
‘tetragamma’, ‘pentagamma’, and ‘hexagamma’.
The functions choose and lchoose return binomial
coefficients and the logarithms of their absolute values. Note that
choose(n, k) is defined for all real numbers n and integer
k. For k ≥ 1 it is defined as
n(n-1)…(n-k+1) / k!,
as 1 for k = 0 and as 0 for negative k.
Non-integer values of k are rounded to an integer, with a warning.
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 internal generic primitive functions: methods can be
defined for them individually or via the
Math group generic.
Source
gamma, lgamma, beta and lbeta are based on
C translations of Fortran subroutines by W. Fullerton of Los Alamos
Scientific Laboratory (now available as part of SLATEC).
digamma, trigamma and psigamma are based on
Amos, D. E. (1983). A portable Fortran subroutine for
derivatives 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. (For gamma and lgamma.)
Abramowitz, M. and Stegun, I. A. (1972)
Handbook of Mathematical Functions. New York: Dover.
https://en.wikipedia.org/wiki/Abramowitz_and_Stegun provides
links to the full text which is in public domain.
Chapter 6: Gamma and Related Functions.
See Also
Arithmetic for simple, sqrt for
miscellaneous mathematical functions and Bessel for the
real Bessel functions.
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.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: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)
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(base)
> png(filename="/home/ddbj/snapshot/RGM3/R_rel/result/base/Special.Rd_%03d_medium.png", width=480, height=480)
> ### Name: Special
> ### Title: Special Functions of Mathematics
> ### Aliases: Special beta lbeta gamma lgamma psigamma digamma trigamma
> ### choose lchoose factorial lfactorial
> ### Keywords: math
>
> ### ** Examples
>
> require(graphics)
>
> choose(5, 2)
[1] 10
> for (n in 0:10) print(choose(n, k = 0:n))
[1] 1
[1] 1 1
[1] 1 2 1
[1] 1 3 3 1
[1] 1 4 6 4 1
[1] 1 5 10 10 5 1
[1] 1 6 15 20 15 6 1
[1] 1 7 21 35 35 21 7 1
[1] 1 8 28 56 70 56 28 8 1
[1] 1 9 36 84 126 126 84 36 9 1
[1] 1 10 45 120 210 252 210 120 45 10 1
>
> factorial(100)
[1] 9.332622e+157
> lfactorial(10000)
[1] 82108.93
>
> ## 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.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")
+ }
-4 : 0 0 1 -4 10 -20
-3 : 0 0 1 -3 6 -10
-2 : 0 0 1 -2 3 -4
-1 : 0 0 1 -1 1 -1
0 : 0 0 1 0 0 0
1 : 0 0 1 1 0 0
2 : 0 0 1 2 1 0 0
3 : 0 0 1 3 3 1 0 0
4 : 0 0 1 4 6 4 1 0 0
5 : 0 0 1 5 10 10 5 1 0 0
6 : 0 0 1 6 15 20 15 6 1 0 0
7 : 0 0 1 7 21 35 35 21 7 1 0 0
8 : 0 0 1 8 28 56 70 56 28 8 1 0 0
9 : 0 0 1 9 36 84 126 126 84 36 9 1 0 0
10 : 0 0 1 10 45 120 210 252 210 120 45 10 1 0 0
>
> ## 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))
k
[1,] -1 0.00000000 0.00000000
[2,] 0 1.00000000 1.00000000
[3,] 1 0.50000000 0.50000000
[4,] 2 -0.12500000 -0.12500000
[5,] 3 0.06250000 0.06250000
[6,] 4 -0.03906250 -0.03906250
[7,] 5 0.02734375 0.02734375
[8,] 6 -0.02050781 -0.02050781
Warning message:
In gamma(x + 1) : NaNs produced
>
> ## 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)
[1] 1.118034
> sum(choose(1/2, k)* .25^k)
[1] 1.118034
>
> ## Don't show:
> k. <- 1:9
> stopifnot(all.equal( (choose(1/2, k.) -> ck.),
+ mychoose(1/2, k.)),
+ all.equal(lchoose(1/2, k.), log(abs(ck.))),
+ all.equal(sqrt(1.25),
+ sum(choose(1/2, k)* .25^k)))
> ## End(Don't show)
>
>
>
>
>
> dev.off()
null device
1
>