与beta函数和gamma函数有关的特殊数学函数。
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 : 非负数值向量。
x, n : 数字向量。
k, deriv : 整数向量。
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)