Special {base} | R Documentation |
Special mathematical functions related to the beta and gamma functions.
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)
a , b |
non-negative numeric vectors. |
x , n |
numeric vectors. |
k , deriv |
integer vectors. |
The functions beta
and lbeta
return the beta function
and the natural logarithm of the beta function,
B(a,b) = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}.
The formal definition is
B(a, b) = \int_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
\Gamma(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)
\Gamma(x) = \int_0^\infty t^{a-1} e^{-t} dt
for all real x
except zero and negative integers (when
NaN
is returned).
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 \psi(x)
.
\code{digamma(x)} = \psi(x) = \frac{d}{dx}\ln\Gamma(x) =
\frac{\Gamma'(x)}{\Gamma(x)}
This is often called the ‘polygamma’ function, e.g. in
Abramowitz and Stegun (section 6.4.1, page 260); and its higher
derivatives (deriv = 2:4
) have occasionally been called
‘tetragamma’, ‘pentagamma’, and ‘hexagamma’.
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
\ge 1
as n(n-1)\cdots(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 generic: methods can be defined for them individually or
via the Math
group generic.
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 Software 9(4), 494–502.
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.
Arithmetic
for simple, sqrt
for
miscellaneous mathematical functions and Bessel
for the
real Bessel functions.
For the incomplete gamma function see pgamma
.
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(.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 <- 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)