Abstract
最近接触到 von Mises–Fisher distribution, 其概率密度如下:
f
p
(
x
;
μ
,
κ
)
=
κ
p
2
−
1
(
2
π
)
p
2
I
p
2
−
1
(
κ
)
e
κ
μ
⊺
x
\begin{aligned} f_{p}(\bm{x}; \bm{\mu}, \kappa) = \frac{\kappa^{\frac{p}{2}-1}} {(2\pi)^{\frac{p}{2}} I_{\frac{p}{2}-1}(\kappa)} e^{\kappa\bm{\mu^\intercal}\bm{x}} \end{aligned}
fp(x;μ,κ)=(2π)2pI2p−1(κ)κ2p−1eκμ⊺x 其中
I
α
(
x
)
I_{\alpha}(x)
Iα(x) 是
α
\alpha
α 阶第一类贝塞尔函数.查阅 Wikipedia, 其公式为:
前者是定义式, 后者是积分式. 其计算比较麻烦, 更不用说求导了, 好在这类带有
e
∗
∗
∗
e^{***}
e∗∗∗ 的积分式能推出递推式, 主要有:
1. 推导 d I α d x = I α − 1 ( x ) − α x I α ( x ) \frac{dI_{\alpha}}{dx} = I_{\alpha-1}(x) - \frac{\alpha}{x} I_{\alpha}(x) dxdIα=Iα−1(x)−xαIα(x):
I α ( x ) = ∑ m = 0 ∞ ( x 2 ) 2 m + α m ! Γ ( m + α + 1 ) d I α d x = ∑ m = 0 ∞ ( 2 m + α ) ( x 2 ) 2 m + α − 1 2 m ! Γ ( m + α + 1 ) = ∑ m = 0 ∞ ( 2 m + 2 α − α ) ( x 2 ) 2 m + α − 1 2 m ! Γ ( m + α + 1 ) = ∑ m = 0 ∞ ( m + α ) ( x 2 ) 2 m + [ α − 1 ] m ! ( m + α ) Γ ( m + [ α − 1 ] + 1 ) − α x ∑ m = 0 ∞ ( x 2 ) ( x 2 ) 2 m + α − 1 m ! Γ ( m + α + 1 ) = I α − 1 ( x ) − α x I α ( x ) \begin{aligned} I_{\alpha}(x) &= \sum_{m=0}^{\infin} \frac{(\frac{x}{2})^{2m+\alpha}}{m!\Gamma(m+\alpha+1)} \\ \frac{dI_{\alpha}}{dx} &= \sum_{m=0}^{\infin} \frac{(2m+\alpha)(\frac{x}{2})^{2m+\alpha-1}}{2m!\Gamma(m+\alpha+1)} \\ &= \sum_{m=0}^{\infin} \frac{(2m+2\alpha-\alpha)(\frac{x}{2})^{2m+\alpha-1}}{2m!\Gamma(m+\alpha+1)} \\ &= \sum_{m=0}^{\infin} \frac{(m+\alpha)(\frac{x}{2})^{2m+[\alpha-1]}}{m!(m+\alpha)\Gamma(m+[\alpha-1]+1)} - \frac{\alpha}{x}\sum_{m=0}^{\infin} \frac{(\frac{x}{2})(\frac{x}{2})^{2m+\alpha-1}}{m!\Gamma(m+\alpha+1)} \\ &= I_{\alpha-1}(x) - \frac{\alpha}{x} I_{\alpha}(x) \end{aligned} Iα(x)dxdIα=m=0∑∞m!Γ(m+α+1)(2x)2m+α=m=0∑∞2m!Γ(m+α+1)(2m+α)(2x)2m+α−1=m=0∑∞2m!Γ(m+α+1)(2m+2α−α)(2x)2m+α−1=m=0∑∞m!(m+α)Γ(m+[α−1]+1)(m+α)(2x)2m+[α−1]−xαm=0∑∞m!Γ(m+α+1)(2x)(2x)2m+α−1=Iα−1(x)−xαIα(x)
2. 推导 d I α d x = I α + 1 ( x ) + α x I α ( x ) \frac{dI_{\alpha}}{dx} = I_{\alpha+1}(x) + \frac{\alpha}{x} I_{\alpha}(x) dxdIα=Iα+1(x)+xαIα(x):
d
I
α
d
x
=
∑
m
=
0
∞
(
2
m
+
α
)
(
x
2
)
2
m
+
α
−
1
2
m
!
Γ
(
m
+
α
+
1
)
=
∑
m
=
0
∞
m
(
x
2
)
2
m
+
α
−
1
m
!
Γ
(
m
+
α
+
1
)
+
∑
m
=
0
∞
α
(
x
2
)
2
m
+
α
−
1
2
m
!
Γ
(
m
+
α
+
1
)
\begin{aligned} \frac{dI_{\alpha}}{dx} &= \sum_{m=0}^{\infin} \frac{(2m+\alpha)(\frac{x}{2})^{2m+\alpha-1}}{2m!\Gamma(m+\alpha+1)} \\ &= \sum_{m=0}^{\infin} \frac{m(\frac{x}{2})^{2m+\alpha-1}}{m!\Gamma(m+\alpha+1)} + \sum_{m=0}^{\infin} \frac{\alpha(\frac{x}{2})^{2m+\alpha-1}}{2m!\Gamma(m+\alpha+1)} \\ \end{aligned}
dxdIα=m=0∑∞2m!Γ(m+α+1)(2m+α)(2x)2m+α−1=m=0∑∞m!Γ(m+α+1)m(2x)2m+α−1+m=0∑∞2m!Γ(m+α+1)α(2x)2m+α−1 后一项就是
α
x
I
α
(
x
)
\frac{\alpha}{x} I_{\alpha}(x)
xαIα(x), 但前面一项看着不好搞, 式中
x
2
\frac{x}{2}
2x 的次数是
(
2
m
+
α
−
1
)
(2m+\alpha-1)
(2m+α−1), 却能推到
I
α
+
1
(
x
)
I_{\alpha+1}(x)
Iα+1(x), 起初还从
Γ
(
m
+
α
+
1
)
=
Γ
(
m
+
[
α
+
1
]
+
1
)
m
+
α
+
1
\Gamma(m+\alpha+1) = \frac{\Gamma(m+[\alpha+1]+1)}{m+\alpha+1}
Γ(m+α+1)=m+α+1Γ(m+[α+1]+1) 搞, 但发现结果是:
∑
m
=
0
∞
m
(
x
2
)
2
m
+
α
−
1
m
!
Γ
(
m
+
α
+
1
)
=
∑
m
=
0
∞
m
(
m
+
α
+
1
)
(
x
2
)
2
m
+
α
−
1
m
!
Γ
(
m
+
[
α
+
1
]
+
1
)
\sum_{m=0}^{\infin} \frac{m(\frac{x}{2})^{2m+\alpha-1}}{m!\Gamma(m+\alpha+1)} = \sum_{m=0}^{\infin} \frac{m(m+\alpha+1)(\frac{x}{2})^{2m+\alpha-1}}{m!\Gamma(m+[\alpha+1]+1)}
m=0∑∞m!Γ(m+α+1)m(2x)2m+α−1=m=0∑∞m!Γ(m+[α+1]+1)m(m+α+1)(2x)2m+α−1 次数都不对, 咋搞? 发现若约去
m
m
m, 则
(
m
−
1
)
!
∣
m
=
0
=
(
−
1
)
!
(m-1)!|_{m=0} = (-1)!
(m−1)!∣m=0=(−1)!, 那就单独拿出来看看:
m
(
x
2
)
2
m
+
α
−
1
m
!
Γ
(
m
+
α
+
1
)
∣
m
=
0
=
0
\frac{m(\frac{x}{2})^{2m+\alpha-1}}{m!\Gamma(m+\alpha+1)}|_{m=0} = 0
m!Γ(m+α+1)m(2x)2m+α−1∣m=0=0, 则:
d
I
α
d
x
=
∑
m
=
0
∞
m
(
x
2
)
2
m
+
α
−
1
m
!
Γ
(
m
+
α
+
1
)
+
α
x
I
α
(
x
)
=
∑
m
=
1
∞
(
x
2
)
2
m
+
α
−
1
(
m
−
1
)
!
Γ
(
m
+
α
+
1
)
+
α
x
I
α
(
x
)
t
=
m
−
1
⇒
=
∑
t
=
0
∞
(
x
2
)
2
(
t
+
1
)
+
α
−
1
t
!
Γ
(
(
t
+
1
)
+
α
+
1
)
+
α
x
I
α
(
x
)
=
∑
t
=
0
∞
(
x
2
)
2
t
+
[
α
+
1
]
t
!
Γ
(
t
+
[
α
+
1
]
+
1
)
+
α
x
I
α
(
x
)
=
I
α
+
1
(
x
)
+
α
x
I
α
(
x
)
\begin{aligned} \frac{dI_{\alpha}}{dx} &= \sum_{m=0}^{\infin} \frac{m(\frac{x}{2})^{2m+\alpha-1}}{m!\Gamma(m+\alpha+1)} + \frac{\alpha}{x} I_{\alpha}(x) \\ &= \sum_{m=1}^{\infin} \frac{(\frac{x}{2})^{2m+\alpha-1}}{(m-1)!\Gamma(m+\alpha+1)} + \frac{\alpha}{x} I_{\alpha}(x) \\ t = m-1 \Rightarrow &= \sum_{t=0}^{\infin} \frac{(\frac{x}{2})^{2(t+1)+\alpha-1}}{t!\Gamma((t+1)+\alpha+1)} + \frac{\alpha}{x} I_{\alpha}(x) \\ &= \sum_{t=0}^{\infin} \frac{(\frac{x}{2})^{2t+[\alpha+1]}}{t!\Gamma(t+[\alpha+1]+1)} + \frac{\alpha}{x} I_{\alpha}(x) \\ &= I_{\alpha+1}(x) + \frac{\alpha}{x} I_{\alpha}(x) \end{aligned}
dxdIαt=m−1⇒=m=0∑∞m!Γ(m+α+1)m(2x)2m+α−1+xαIα(x)=m=1∑∞(m−1)!Γ(m+α+1)(2x)2m+α−1+xαIα(x)=t=0∑∞t!Γ((t+1)+α+1)(2x)2(t+1)+α−1+xαIα(x)=t=0∑∞t!Γ(t+[α+1]+1)(2x)2t+[α+1]+xαIα(x)=Iα+1(x)+xαIα(x)
3. 其他两个递推式
推了这两个, 那么两式相加得: I α − 1 ( x ) + I α + 1 ( x ) = 2 d I α d x I_{\alpha-1}(x) + I_{\alpha+1}(x) = 2 \frac{dI_{\alpha}}{dx} Iα−1(x)+Iα+1(x)=2dxdIα 两式相减得: I α − 1 ( x ) − I α + 1 ( x ) = 2 I α I_{\alpha-1}(x) - I_{\alpha+1}(x) = 2 I_{\alpha} Iα−1(x)−Iα+1(x)=2Iα 就很自然了.
4. 积分式的递推式证明
见第一类修正贝塞尔函数的递推公式怎么证明?, 这是阶数
α
\alpha
α 为整数时的情况.
5. Exponentially Scaled Modified Bessel Function of the First Kind.
Defined as::
ive(v, z) = iv(v, z) * exp(-abs(z.real))
所以当假定 z . r e a l ≥ 0 z.real \ge 0 z.real≥0 时, 求导数就按照莱布尼茨法则, 就可以了:
ive(v, z)' = ive(ctx.v - 1, z) - ive(ctx.v, z) * (ctx.v + z) / z
6. 用 scipy.special
包进行验证
import numpy as np
from scipy import special
v = 5.5
z = 3
iv = special.iv(v, z)
ive = special.ive(v, z) # iv * np.exp(-abs(Re(z)))
print(iv) # 0.04532335799965591
print(ive) # 0.0022565171233900425
print(iv * np.exp(-z)) # 验证 ive = iv * exp(-z), 0.002256517123390042
# %% 验证导数
ivp = special.ivp(v, z)
ivp_1 = (special.iv(v - 1, z) + special.iv(v + 1, z)) / 2
ivp_2 = special.iv(v - 1, z) - v / z * iv
ivp_3 = special.iv(v + 1, z) + v / z * iv
print(ivp_1) # 0.09310529508597687
print(ivp_2) # 0.09310529508597686
print(ivp_3) # 0.09310529508597688
# %% 验证ive的导数
ivep = (ivp - iv) * np.exp(-z)
ivep_1 = special.ive(v - 1, z) - (z + v) / z * ive
print(ivep) # 0.0023789225684656356
print(ivep_1) # 0.002378922568465644