34.28 特殊函数

34.28.1 阶乘

  • 阶乘 \(n! = 1\times 2\times 3\cdots n\)
  • 双阶乘 \((2n+1)!! = 1 \times 3\times 5 \times \cdots \times (2n+1), n = 0,1,2,\cdots\)
factorial(5) # 阶乘
## [1] 120
seq(from = 1, to = 5, length.out = 3)
## [1] 1 3 5
prod(seq(from = 1, to = 5, length.out = 3)) # 连乘 双阶乘
## [1] 15
seq(5)
## [1] 1 2 3 4 5
cumprod(seq(5)) # 累积
## [1]   1   2   6  24 120
cumsum(seq(5)) # 累和
## [1]  1  3  6 10 15

此外还有 cummaxcummin

  • 组合数 \(C_{n}^{k} = \frac{n(n-1)…(n-k+1)}{k!}\)

\(C_{5}^{3} = \frac{5 \times 4 \times 3}{3 \times 2 \times 1}\)

choose(5,3)
## [1] 10

斯特林公式

34.28.2 伽马函数

\(\Gamma(x) = \int_{0}^{\infty} t^{x-1}\exp(-t)dt\) \(\Gamma(n) = (n-1)!, n \in \mathbb{Z}^{+}\)

gamma(2)
## [1] 1
gamma(10)
## [1] 362880
gamma2 <- function(t,x){
  t^(x-1)*exp(-t)
}
integrate(gamma2, lower = 0, upper = + Inf, x = 10)
## 362880 with absolute error < 0.025
  • psigamma(x, deriv) 表示 \(\psi(x)\)deriv 阶导数

\(\mathrm{digamma}(x) \triangleq \psi(x) = \frac{d}{dx}{\ln \Gamma(x)} = \Gamma'(x) / \Gamma(x)\)

# 例1
x <- 2
eval(deriv(~ gamma(x), "x"))/gamma(x)
## [1] 1
## attr(,"gradient")
##              x
## [1,] 0.4227843
# 与此等价
psigamma(2, 0)
## [1] 0.4227843
digamma(x) # psi(x) 的一阶导数
## [1] 0.4227843
trigamma(x) # psi(x) 的二阶导数
## [1] 0.6449341
# 例2
eval(deriv(~ psigamma(x, 1), "x"))
## [1] 0.6449341
## attr(,"gradient")
##               x
## [1,] -0.4041138
# 与此等价
psigamma(2, 2)
## [1] -0.4041138
# 注意与下面这个例子比较
dx2x <- deriv(~ x^3, "x")
eval(dx2x)
## [1] 8
## attr(,"gradient")
##       x
## [1,] 12

34.28.3 贝塔函数

\(B(a,b) = \Gamma(a)\Gamma(b)/\Gamma(a+b) = \int_{0}^{1} t^{a-1} (1-t)^{b-1} dt\)

beta(1, 1)
## [1] 1
beta(2, 3)
## [1] 0.08333333
beta2 <- function(t, a, b) {
  t^(a - 1) * (1 - t)^(b - 1)
}
integrate(beta2, lower = 0, upper = 1, a = 2, b = 3)
## 0.08333333 with absolute error < 9.3e-16

34.28.4 贝塞尔函数

besselI(x, nu, expon.scaled = FALSE) # 修正的第一类
besselK(x, nu, expon.scaled = FALSE) # 修正的第二类
besselJ(x, nu) # 第一类
besselY(x, nu) # 第二类
  • \(\nu\) 贝塞尔函数的阶,可以是分数
  • expon.scaled 是否使用指数表示
nus <- c(0:5, 10, 20)
x <- seq(0, 4, length.out = 501)
plot(x, x,
  ylim = c(0, 6), ylab = "", type = "n",
  main = "Bessel Functions  I_nu(x)"
)
for (nu in nus) lines(x, besselI(x, nu = nu), col = nu + 2)
legend(0, 6, legend = paste("nu=", nus), col = nus + 2, lwd = 1)