Modified Bessel Function of the First Kind

Abstract

最近接触到 von Mises–Fisher distribution, 其概率密度如下:
$$\begin{array}{r}{f}_p(\mathbf{x};\mathbf{\mu },\kappa )=\frac{{\kappa }^{\frac{p}{2}-1}}{(2\pi {)}^{\frac{p}{2}}{I}_{\frac{p}{2}-1}(\kappa )}{e}^{\kappa {\mathbf{\mu }}^{⊺}\mathbf{x}}\end{array}$$ 其中 $I_{\alpha }(x)$ 是 $\alpha$ 阶第一类贝塞尔函数.查阅 Wikipedia, 其公式为:


前者是定义式, 后者是积分式. 其计算比较麻烦, 更不用说求导了, 好在这类带有 $e^{\ast \ast \ast }$ 的积分式能推出递推式, 主要有:

1. 推导 $\frac{d{I}_{\alpha }}{dx}=I_{\alpha -1}(x)-\frac{\alpha }{x}{I}_{\alpha }(x)$:

$$\begin{array}{rl}{I}_{\alpha }(x)& =\sum _{m=0}^{\infty }\frac{(\frac{x}{2}{)}^{2m+\alpha }}{m!\Gamma (m+\alpha +1)}\\ \frac{d{I}_{\alpha }}{dx}& =\sum _{m=0}^{\infty }\frac{(2m+\alpha )(\frac{x}{2}{)}^{2m+\alpha -1}}{2m!\Gamma (m+\alpha +1)}\\ & =\sum _{m=0}^{\infty }\frac{(2m+2\alpha -\alpha )(\frac{x}{2}{)}^{2m+\alpha -1}}{2m!\Gamma (m+\alpha +1)}\\ & =\sum _{m=0}^{\infty }\frac{(m+\alpha )(\frac{x}{2}{)}^{2m+[\alpha -1]}}{m!(m+\alpha )\Gamma (m+[\alpha -1]+1)}-\frac{\alpha }{x}\sum _{m=0}^{\infty }\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{array}$$

2. 推导 $\frac{d{I}_{\alpha }}{dx}=I_{\alpha +1}(x)+\frac{\alpha }{x}{I}_{\alpha }(x)$:

$$\begin{array}{rl}\frac{d{I}_{\alpha }}{dx}& =\sum _{m=0}^{\infty }\frac{(2m+\alpha )(\frac{x}{2}{)}^{2m+\alpha -1}}{2m!\Gamma (m+\alpha +1)}\\ & =\sum _{m=0}^{\infty }\frac{m(\frac{x}{2}{)}^{2m+\alpha -1}}{m!\Gamma (m+\alpha +1)}+\sum _{m=0}^{\infty }\frac{\alpha (\frac{x}{2}{)}^{2m+\alpha -1}}{2m!\Gamma (m+\alpha +1)}\end{array}$$ 后一项就是 $\frac{\alpha }{x}{I}_{\alpha }(x)$, 但前面一项看着不好搞, 式中 $\frac{x}{2}$ 的次数是 $(2m+\alpha -1)$, 却能推到 $I_{\alpha +1}(x)$, 起初还从 $\Gamma (m+\alpha +1)=\frac{\Gamma (m+[\alpha +1]+1)}{m+\alpha +1}$ 搞, 但发现结果是: $$\sum _{m=0}^{\infty }\frac{m(\frac{x}{2}{)}^{2m+\alpha -1}}{m!\Gamma (m+\alpha +1)}=\sum _{m=0}^{\infty }\frac{m(m+\alpha +1)(\frac{x}{2}{)}^{2m+\alpha -1}}{m!\Gamma (m+[\alpha +1]+1)}$$ 次数都不对, 咋搞? 发现若约去 $m$, 则 $(m-1)!{|}_{m=0}=(-1)!$, 那就单独拿出来看看: $\frac{m(\frac{x}{2}{)}^{2m+\alpha -1}}{m!\Gamma (m+\alpha +1)}{|}_{m=0}=0$, 则:
$$\begin{array}{rl}\frac{d{I}_{\alpha }}{dx}& =\sum _{m=0}^{\infty }\frac{m(\frac{x}{2}{)}^{2m+\alpha -1}}{m!\Gamma (m+\alpha +1)}+\frac{\alpha }{x}{I}_{\alpha }(x)\\ & =\sum _{m=1}^{\infty }\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}^{\infty }\frac{(\frac{x}{2}{)}^{2(t+1)+\alpha -1}}{t!\Gamma ((t+1)+\alpha +1)}+\frac{\alpha }{x}{I}_{\alpha }(x)\\ & =\sum _{t=0}^{\infty }\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{array}$$

3. 其他两个递推式

推了这两个, 那么两式相加得: $$I_{\alpha -1}(x)+I_{\alpha +1}(x)=2\frac{d{I}_{\alpha }}{dx}$$ 两式相减得: $$I_{\alpha -1}(x)-I_{\alpha +1}(x)=2I_{\alpha }$$ 就很自然了.

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.real\ge 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

补充材料