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)=\frac{2\alpha }{x}{I}_{\alpha }$$ 就很自然了.

4. 积分式的递推式证明

见第一类修正贝塞尔函数的递推公式怎么证明？, 这是阶数 $\alpha$ 为整数时($sin(\alpha \pi )=0$)的情况.

那 $\alpha$ 不为整数时呢? 发现 $$\begin{array}{rl}-sin(nx)e^{\alpha cosx}{|}_0^{\pi }=-sin(n\pi )e^{-\alpha }& \ne 0\\ \frac{sin(\alpha \pi )}{\pi }& \ne 0\end{array}$$ 那么我们就来推一遍。由于太长，我们分前后两部分。

4.1 递推式 ①

前半部分：$$\begin{array}{rl}& \frac{1}{\pi }\left({\int }_0^{\pi }{e}^{xcos\theta }cos((\alpha -1)\theta )d\theta -{\int }_0^{\pi }{e}^{xcos\theta }cos((\alpha +1)\theta )d\theta \right)\\ =& \frac{1}{\pi }{\int }_0^{\pi }{e}^{xcos\theta }\left(cos((\alpha -1)\theta )-cos((\alpha +1)\theta )\right)d\theta \\ =& \frac{2}{\pi }{\int }_0^{\pi }sin\theta sin(\alpha \theta )e^{xcos\theta }d\theta \\ =& -\frac{2}{\pi x}{\int }_0^{\pi }sin(\alpha \theta )d(e^{xcos\theta })\\ =& -\frac{2}{\pi x}\left(sin(\alpha \theta )e^{xcos\theta }{|}_0^{\pi }-{\int }_0^{\pi }{e}^{xcos\theta }d(sin(\alpha \theta ))\right)\\ =& -\frac{2sin(\alpha \theta )}{\pi x}{e}^{xcos\theta }{|}_0^{\pi }+\frac{2\alpha }{x}\frac{1}{\pi }{\int }_0^{\pi }{e}^{xcos\theta }cos(\alpha \theta )d\theta \end{array}$$ 再看积分式的后半部分: $$\begin{array}{rl}& \frac{sin((\alpha -1)\pi )}{\pi }{\int }_0^{\infty }{e}^{-xcosht-(\alpha -1)t}dt-\frac{sin((\alpha +1)\pi )}{\pi }{\int }_0^{\infty }{e}^{-xcosht-(\alpha +1)t}dt\\ =& -\frac{sin(\alpha \pi )}{\pi }{\int }_0^{\infty }\left(e^{-xcosht-(\alpha -1)t}-e^{-xcosht-(\alpha +1)t}\right)dt\\ =& -\frac{sin(\alpha \pi )}{\pi }{\int }_0^{\infty }{e}^{-xcosht-\alpha t}\left(e^t-e^{-t}\right)dt\\ =& -2\frac{sin(\alpha \pi )}{\pi }{\int }_0^{\infty }sinht{e}^{-xcosht-\alpha t}dt\\ =& \frac{2}{x}\frac{sin(\alpha \pi )}{\pi }{\int }_0^{\infty }{e}^{-\alpha t}d\left(e^{-xcosht}\right)\\ =& \frac{2}{x}\frac{sin(\alpha \pi )}{\pi }\left(e^{-xcosht-\alpha t}{|}_0^{\infty }+{\int }_0^{\infty }\alpha e^{-xcosht-\alpha t}dt\right)\\ =& \frac{2sin(\alpha \pi )}{\pi x}{e}^{-xcosht-\alpha t}{|}_0^{\infty }+\frac{2\alpha }{x}\frac{sin(\alpha \pi )}{\pi }{\int }_0^{\infty }{e}^{-xcosht-\alpha t}dt\end{array}$$ 此时, 就看 $-\frac{2sin(\alpha \theta )}{\pi x}{e}^{xcos\theta }{|}_0^{\pi }-\frac{2sin(\alpha \pi )}{\pi x}{e}^{-xcosht-\alpha t}{|}_0^{\infty }$ 了：$$\begin{array}{rl}& -\frac{2sin(\alpha \theta )}{\pi x}{e}^{xcos\theta }{|}_0^{\pi }-\frac{2sin(\alpha \pi )}{\pi x}{e}^{-xcosht-\alpha t}{|}_0^{\infty }\\ =& -\frac{2sin(\alpha \pi )}{\pi x}{e}^{-x}-\frac{2sin(\alpha \pi )}{\pi x}(-e^{-x})\\ =& -\frac{2sin(\alpha \pi )}{\pi x}{e}^{-x}+\frac{2sin(\alpha \pi )}{\pi x}{e}^{-x}\\ =& 0\end{array}$$ $I_{\alpha -1}(x)-I_{\alpha +1}(x)=\frac{2\alpha }{x}{I}_{\alpha }$ 得证。

4.2 递推式 ②

前半部分：$$\begin{array}{rl}& \frac{1}{\pi }\left({\int }_0^{\pi }{e}^{xcos\theta }cos((\alpha -1)\theta )d\theta +{\int }_0^{\pi }{e}^{xcos\theta }cos((\alpha +1)\theta )d\theta \right)\\ =& \frac{1}{\pi }{\int }_0^{\pi }{e}^{xcos\theta }\left(cos((\alpha -1)\theta )+cos((\alpha +1)\theta )\right)d\theta \\ =& \frac{2}{\pi }{\int }_0^{\pi }cos\theta cos(\alpha \theta )e^{xcos\theta }d\theta \\ =& \frac{2}{\pi }\frac{d\left({\int }_0^{\pi }cos(\alpha \theta )e^{xcos\theta }d\theta \right)}{dx}\end{array}$$ 再看积分式的后半部分: $$\begin{array}{rl}& \frac{sin((\alpha -1)\pi )}{\pi }{\int }_0^{\infty }{e}^{-xcosht-(\alpha -1)t}dt+\frac{sin((\alpha +1)\pi )}{\pi }{\int }_0^{\infty }{e}^{-xcosht-(\alpha +1)t}dt\\ =& -\frac{sin(\alpha \pi )}{\pi }{\int }_0^{\infty }\left(e^{-xcosht-(\alpha -1)t}+e^{-xcosht-(\alpha +1)t}\right)dt\\ =& -\frac{sin(\alpha \pi )}{\pi }{\int }_0^{\infty }{e}^{-xcosht-\alpha t}\left(e^t+e^{-t}\right)dt\\ =& 2\frac{sin(\alpha \pi )}{\pi }{\int }_0^{\infty }(-cosht)e^{-xcosht-\alpha t}dt\\ =& 2\frac{sin(\alpha \pi )}{\pi }\frac{d({\int }_0^{\infty }{e}^{-xcosht-\alpha t}dt)}{dx}\end{array}$$ 合并: $$\begin{array}{rl}前半部分-后半部分=& \frac{2}{\pi }\frac{d\left({\int }_0^{\pi }cos(\alpha \theta )e^{xcos\theta }d\theta \right)}{dx}-2\frac{sin(\alpha \pi )}{\pi }\frac{d({\int }_0^{\infty }{e}^{-xcosht-\alpha t}dt)}{dx}\\ =& 2\left(\frac{d(\frac{1}{\pi }{\int }_0^{\pi }cos(\alpha \theta )e^{xcos\theta }-\frac{sin(\alpha \pi )}{\pi }{\int }_0^{\infty }{e}^{-xcosht-\alpha t}dt)}{dx}\right)\\ =& 2\frac{d{I}_{\alpha }(x)}{dx}\end{array}$$ $I_{\alpha -1}(x)+I_{\alpha +1}(x)=2\frac{d{I}_{\alpha }(x)}{dx}$ 得证。

其他两个递推式可由上面两个式子推导出来: ①+② 可得 ③, ①-② 可得 ④.

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

至于为什么要用 exp(-abs(z.real)) 来缩放 iv 函数, 大概是因为 iv 这个函数是呈指数增长的吧:

Wikipedia: Modified Bessel Functions of the First Kind, Iα(x), for α = 0, 1, 2, 3

当经过 $exp(-x)$ 的缩放后, 函数应该会变得更平滑吧.

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

7. 补充材料

7.1 公式之间的关系

$$\begin{array}{rl}{I}_{\alpha }(x)=& \frac{(\frac{x}{2}{)}^{\alpha }}{\sqrt{\pi }\Gamma (\alpha +\frac{1}{2})}{\int }_{-1}^1{e}^{-xt}(1-t^2{)}^{\alpha -\frac{1}{2}}dt(Re\alpha >-\frac{1}{2})\\ =& \frac{(\frac{x}{2}{)}^{\alpha }}{\sqrt{\pi }\Gamma (\alpha +\frac{1}{2})}{\int }_0^{\pi }cosh(xcos\theta )si{n}^{2\alpha }\theta d\theta (Re\alpha >-\frac{1}{2})\end{array}$$ 如何转换? 首先想到的是 $t=-cos\theta$: $$\begin{array}{rl}{I}_{\alpha }(x)=& \frac{(\frac{x}{2}{)}^{\alpha }}{\sqrt{\pi }\Gamma (\alpha +\frac{1}{2})}{\int }_0^{\pi }{e}^{xcos\theta }si{n}^{2\alpha -1}\theta (sin\theta d\theta )\\ =& \frac{(\frac{x}{2}{)}^{\alpha }}{\sqrt{\pi }\Gamma (\alpha +\frac{1}{2})}{\int }_0^{\pi }{e}^{xcos\theta }si{n}^{2\alpha }\theta d\theta \end{array}\text{\hspace{1em}(1)}$$ 再令 $t=cos\theta$: $$\begin{array}{rl}{I}_{\alpha }(x)=& \frac{(\frac{x}{2}{)}^{\alpha }}{\sqrt{\pi }\Gamma (\alpha +\frac{1}{2})}{\int }_{\pi }^0{e}^{-xcos\theta }si{n}^{2\alpha -1}\theta (-sin\theta d\theta )\\ =& \frac{(\frac{x}{2}{)}^{\alpha }}{\sqrt{\pi }\Gamma (\alpha +\frac{1}{2})}{\int }_0^{\pi }{e}^{-xcos\theta }si{n}^{2\alpha }\theta d\theta \end{array}\text{\hspace{1em}(2)}$$ 两种变换求平均, 就得到了.