von Mises-Fisher Distribution

1. 概率密度函数

关于 Bessel Function 请参阅《贝塞尔函数-Wikipedia》《Modified Bessel Function of the First Kind》。

小结: 集中参数 $\kappa \ge 0$ 越大, 方向的分布越集中于平均方向 $\mathbf{\mu }$. 极限地, 当 $\kappa =0$ 时, 方向地均匀分布; 当 $\kappa \to \infty$ 时, 变成 $\mathbf{\mu }$ 处的狄拉克分布.

2. Relation to Normal Distribution

即，将协方差矩阵为 ${\kappa }^{-1}\mathbf{I}$、均值为 $\mathbf{\mu }(\Vert \mathbf{\mu }\Vert =r>0)$ 的 $p$ 元高斯分布的 support 限制在超球 ($\Vert \mathbf{x}\Vert =1$) 上，则变为 von Mises-Fisher Distribution (vMF)。

进一步的推导：$$\begin{array}{rl}{G}_p(\mathbf{x};\mathbf{\mu },\kappa )& ={\left(\sqrt{\frac{\kappa }{2\pi }}\right)}^{p}exp\left(-\kappa \frac{(\mathbf{x}-\mathbf{\mu }{)}^{⊺}(\mathbf{x}-\mathbf{\mu })}{2}\right)\\ & ={\left(\sqrt{\frac{\kappa }{2\pi }}\right)}^{p}exp\left(-\kappa \frac{{\mathbf{x}}^{⊺}\mathbf{x}+{\mathbf{\mu }}^{⊺}\mathbf{\mu }-2{\mathbf{\mu }}^{⊺}\mathbf{x}}{2}\right)\\ & ={\left(\sqrt{\frac{\kappa }{2\pi }}\right)}^{p}exp\left(-\kappa \frac{1+r^2-2{\mathbf{\mu }}^{⊺}\mathbf{x}}{2}\right)\\ & ={\left(\sqrt{\frac{\kappa }{2\pi }}\right)}^{p}exp\left(\frac{-\kappa (1+r^2)}{2}+\kappa {\mathbf{\mu }}^{⊺}\mathbf{x}\right)\\ & =\left[{\left(\sqrt{\frac{\kappa }{2\pi }}\right)}^{p}exp\left(\frac{-\kappa (1+r^2)}{2}\right)\right]exp\left(\kappa r\frac{{\mathbf{\mu }}^{⊺}}{r}\mathbf{x}\right)\end{array}$$ 【注】限制 $\Vert \mathbf{x}\Vert =1$ 后，上式的归一化系数要重新计算：在超球上进行积分。得到的结果为 $f_p(\mathbf{x};r^{-1}\mu ;r\kappa )$。

上图是可视化的二维平面上, 以 $\mathbf{\mu }=(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})$ 为均值、以 $\Sigma =\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$ 为标准差的正态分布, 如果限制分布支撑集为单位圆, 再重新计算归一化系数, 则成了 von-Mises 分布(vMF的二维情形).

小结：von Mises-Fisher Distribution 本质是超球上的正态分布，且各向同性 (isotropic)，协方差矩阵为 ${\kappa }^{-1}\mathbf{I}$。

疑问：有没有不各向同性的 vMF？

3. Estimating the Parameters of the vMF Distribution

根据 Wikipedia, 我们先把 vMF 的概率密度函数表示为: $$\begin{array}{rl}{f}_p(\mathbf{x};\mathbf{\mu },\kappa )=& C_p(\kappa )exp(\kappa {\mathbf{\mu }}^{⊺}\mathbf{x})\\ C_p(\kappa )=& \frac{{\kappa }^{\frac{p}{2}-1}}{(2\pi {)}^{\frac{p}{2}}{I}_{\frac{p}{2}-1}(\kappa )}\end{array}$$ 给定一堆样本 { x 1 , x 2 , ⋯   , x n } \{\bm{x}_1, \bm{x}_2, \cdots, \bm{x}_n\} {x1​,x2​,⋯,xn​}, 独立地采自某 vMF, 现要对其参数 $(\kappa ,\mathbf{\mu })$ 进行最大似然估计: $$\begin{array}{rl}L(\mathbf{\mu };p,\kappa ,\lambda )=& nlog(C_p(\kappa ))+\kappa {\mathbf{\mu }}^{⊺}\sum _{i=1}^n{\mathbf{x}}_i+\lambda (\Vert \mathbf{\mu }\Vert -1)\end{array}$$ 求偏导: $$\begin{array}{rlr}\frac{\partial L}{\partial \mathbf{\mu }}=& \kappa \sum _{i=1}^n{\mathbf{x}}_i+\lambda \frac{\mathbf{\mu }}{\Vert \mathbf{\mu }\Vert }& =0\\ \frac{\partial L}{\partial \kappa }=& n\frac{C_p^{{}^{\prime }}(\kappa )}{C_p(\kappa )}+{\mathbf{\mu }}^{⊺}\sum _{i=1}^n{\mathbf{x}}_i& =0\\ \frac{\partial L}{\partial \lambda }=& \Vert \mathbf{\mu }\Vert -1& =0\end{array}$$ 先看其中比较特别的 $\frac{C_p^{{}^{\prime }}(\kappa )}{C_p(\kappa )}$: $$\begin{array}{rl}{C}_p(\kappa )& =\frac{{\kappa }^{\frac{p}{2}-1}}{(2\pi {)}^{\frac{p}{2}}{I}_{\frac{p}{2}-1}(\kappa )}=(2\pi {)}^{-\frac{p}{2}}\frac{{\kappa }^{\frac{p}{2}-1}}{I_{\frac{p}{2}-1}(\kappa )}\\ C_p^{{}^{\prime }}(\kappa )& =(2\pi {)}^{-\frac{p}{2}}\frac{(\frac{p}{2}-1){\kappa }^{\frac{p}{2}-2}{I}_{\frac{p}{2}-1}(\kappa )-{\kappa }^{\frac{p}{2}-1}{I}_{\frac{p}{2}-1}^{{}^{\prime }}(\kappa )}{I_{\frac{p}{2}-1}^2(\kappa )}\\ \frac{C_p^{{}^{\prime }}(\kappa )}{C_p(\kappa )}& =\frac{(2\pi {)}^{-\frac{p}{2}}\frac{(\frac{p}{2}-1){\kappa }^{\frac{p}{2}-2}{I}_{\frac{p}{2}-1}(\kappa )-{\kappa }^{\frac{p}{2}-1}{I}_{\frac{p}{2}-1}^{{}^{\prime }}(\kappa )}{I_{\frac{p}{2}-1}^2(\kappa )}}{(2\pi {)}^{-\frac{p}{2}}\frac{{\kappa }^{\frac{p}{2}-1}}{I_{\frac{p}{2}-1}(\kappa )}}\\ & =\frac{(\frac{p}{2}-1){\kappa }^{\frac{p}{2}-2}{I}_{\frac{p}{2}-1}(\kappa )-{\kappa }^{\frac{p}{2}-1}{I}_{\frac{p}{2}-1}^{{}^{\prime }}(\kappa )}{{\kappa }^{\frac{p}{2}-1}{I}_{\frac{p}{2}-1}(\kappa )}\\ & =\frac{(\frac{p}{2}-1){\kappa }^{-1}{I}_{\frac{p}{2}-1}(\kappa )-I_{\frac{p}{2}-1}^{{}^{\prime }}(\kappa )}{I_{\frac{p}{2}-1}(\kappa )}\\ & =\frac{(\frac{p}{2}-1){\kappa }^{-1}{I}_{\frac{p}{2}-1}(\kappa )-[I_{\frac{p}{2}}(\kappa )+(\frac{p}{2}-1){\kappa }^{-1}{I}_{\frac{p}{2}-1}(\kappa )]}{I_{\frac{p}{2}-1}(\kappa )}\\ & =-\frac{I_{\frac{p}{2}}(\kappa )}{I_{\frac{p}{2}-1}(\kappa )}<0\end{array}$$ 记 $A_p(\kappa )=-\frac{C_p^{{}^{\prime }}(\kappa )}{C_p(\kappa )}=\frac{I_{\frac{p}{2}}(\kappa )}{I_{\frac{p}{2}-1}(\kappa )}$. 想办法消去前两式中的 $\mathbf{\mu }$: $$\begin{array}{rlr}(\kappa \sum _{i=1}^n{\mathbf{x}}_i{)}^{⊺}(\kappa \sum _{i=1}^n{\mathbf{x}}_i)& =(-\lambda \mathbf{\mu }{)}^{⊺}(-\lambda \mathbf{\mu })& \\ {\kappa }^2(\sum _{i=1}^n{\mathbf{x}}_i{)}^{⊺}(\sum _{i=1}^n{\mathbf{x}}_i)& ={\lambda }^2& \\ \Vert \sum _{i=1}^n{\mathbf{x}}_i\Vert & =\frac{\left|\lambda \right|}{\kappa }& [\lambda 正负未知]\end{array}\text{\hspace{1em}(1)}$$ $$\begin{array}{rlr}\kappa {\mathbf{\mu }}^{⊺}\sum _{i=1}^n{\mathbf{x}}_i+\lambda {\mathbf{\mu }}^{⊺}\mathbf{\mu }& =\kappa {\mathbf{\mu }}^{⊺}\sum _{i=1}^n{\mathbf{x}}_i+\lambda =0& \\ n\kappa \frac{C_p^{{}^{\prime }}(\kappa )}{C_p(\kappa )}+\kappa {\mathbf{\mu }}^{⊺}\sum _{i=1}^n{\mathbf{x}}_i& =\kappa {\mathbf{\mu }}^{⊺}\sum _{i=1}^n{\mathbf{x}}_i+n\kappa \frac{C_p^{{}^{\prime }}(\kappa )}{C_p(\kappa )}=0& \\ n\kappa \frac{C_p^{{}^{\prime }}(\kappa )}{C_p(\kappa )}& =\lambda <0& [\because \frac{C_p^{{}^{\prime }}(\kappa )}{C_p(\kappa )}<0]\end{array}\text{\hspace{1em}(2)}$$ 现在可以去掉 $(1)$ 中的绝对值符号了: $$\begin{array}{rl}\Vert \sum _{i=1}^n{\mathbf{x}}_i\Vert & =-\frac{\lambda }{\kappa }=-n\frac{C_p^{{}^{\prime }}(\kappa )}{C_p(\kappa )}=n{A}_p(\kappa )\\ ⟹A_p(\kappa )& =\frac{1}{n}\Vert \sum _{i=1}^n{\mathbf{x}}_i\Vert \\ & =\Vert \bar{\mathbf{x}}\Vert \end{array}$$ 然而, 关于 $\kappa$ 的方程 $A_p(\kappa )=\Vert \bar{\mathbf{x}}\Vert$ 的求解是相当麻烦的, 一些近似求解方法请参考 Wikipedia. 假装已经求出了 $\kappa$, 然后可以求 $\lambda$: $$\begin{array}{rl}\lambda & =-n\kappa A_p(\kappa )=-n\kappa \Vert \bar{\mathbf{x}}\Vert \end{array}$$ 继续可以求 $\mathbf{\mu }$: $$\begin{array}{rlr}& \kappa \sum _{i=1}^n{\mathbf{x}}_i+\lambda \frac{\mathbf{\mu }}{\Vert \mathbf{\mu }\Vert }=0& \\ ⟶\mathbf{\mu }=& -\frac{\kappa }{\lambda }\sum _{i=1}^n{\mathbf{x}}_i& \\ =& \frac{1}{\Vert {\sum }_{i=1}^n{\mathbf{x}}_i\Vert }\sum _{i=1}^n{\mathbf{x}}_i& [\mathbf{\mu }的估计]\\ =& \frac{1}{n\Vert \bar{\mathbf{x}}\Vert }\sum _{i=1}^n{\mathbf{x}}_i& \\ =& \frac{1}{A_p(\kappa )}\frac{1}{n}\sum _{i=1}^n{\mathbf{x}}_i& \\ ⟶mean=\frac{1}{n}\sum _{i=1}^n{\mathbf{x}}_i=& A_p(\kappa )\mathbf{\mu }& \end{array}$$ 即, 均值(注意不是平均方向)是: $mean=\frac{1}{n}{\sum }_{i=1}^n{\mathbf{x}}_i=A_p(\kappa )\mathbf{\mu }$.

小结: vMF 参数的最大似然估计是: $\{\begin{array}{l}\mathbf{\mu }=\frac{\bar{\mathbf{x}}}{\Vert \bar{\mathbf{x}}\Vert }\\ \kappa =A_p^{-1}(\Vert \bar{\mathbf{x}}\Vert )\end{array}$, 期望值(expected value) 为 $A_p(\kappa )\mathbf{\mu }$.

4. Entropy and KL Divergence

分布 $vMF(\mathbf{\mu },\kappa )$ 的熵为:

$$\begin{array}{rl}⟨-log{f}_p(\mathbf{x};\mathbf{\mu },\kappa ){⟩}_{\mathbf{x}\sim vMF(\mathbf{\mu },\kappa )}& =-{\int }_{球}{f}_p(\mathbf{x};\mathbf{\mu },\kappa )log\left(C_p(\kappa )exp(\kappa {\mathbf{\mu }}^{⊺}\mathbf{x})\right)d\mathbf{x}\\ & =-{\int }_{球}{f}_p(\mathbf{x};\mathbf{\mu },\kappa )\left(log{C}_p(\kappa )+\kappa {\mathbf{\mu }}^{⊺}\mathbf{x}\right)d\mathbf{x}\\ & =-log{C}_p(\kappa ){\int }_{球}{f}_p(\mathbf{x};\mathbf{\mu },\kappa )d\mathbf{x}-\kappa {\mathbf{\mu }}^{⊺}{\int }_{球}\mathbf{x}{f}_p(\mathbf{x};\mathbf{\mu },\kappa )d\mathbf{x}\\ & =-log{C}_p(\kappa )-\kappa {\mathbf{\mu }}^{⊺}{A}_p(\kappa )\mathbf{\mu }\because mean=A_p(\kappa )\mathbf{\mu }\\ & =-log{f}_p(A_p(\kappa )\mathbf{\mu },\mathbf{\mu },\kappa )\\ & =-log{C}_p(\kappa )-\kappa A_p(\kappa )\end{array}$$ 过程中使用了 $$mean={\int }_{球}\mathbf{x}{f}_p(\mathbf{x};\mathbf{\mu },\kappa )d\mathbf{x}=A_p(\kappa )\mathbf{\mu }$$ 是前面讨论参数的最大似然估计时的副产物.

分布 $vMF({\mathbf{\mu }}_0,{\kappa }_0)$ 和 分布 $vMF({\mathbf{\mu }}_1,{\kappa }_1)$ 之间的 KL 散度为:

$$\begin{array}{rl}{⟨log\frac{f_p(\mathbf{x};{\mathbf{\mu }}_0,{\kappa }_0)}{f_p(\mathbf{x};\mathbf{\mu },{\kappa }_1)}⟩}_{\mathbf{x}\sim vMF({\mathbf{\mu }}_0,{\kappa }_0)}& ={\int }_{球}{f}_p(\mathbf{x};{\mathbf{\mu }}_0,{\kappa }_0)log\frac{C_p({\kappa }_0)exp({\kappa }_0{\mathbf{\mu }}_0^{⊺}\mathbf{x})}{C_p({\kappa }_1)exp({\kappa }_1{\mathbf{\mu }}_1^{⊺}\mathbf{x})}d\mathbf{x}\\ & ={\int }_{球}{f}_p(\mathbf{x};{\mathbf{\mu }}_0,{\kappa }_0)(log{C}_p({\kappa }_0)+{\kappa }_0{\mathbf{\mu }}_0^{⊺}\mathbf{x})d\mathbf{x}-{\int }_{球}{f}_p(\mathbf{x};{\mathbf{\mu }}_0,{\kappa }_0)(log{C}_p({\kappa }_1)+{\kappa }_1{\mathbf{\mu }}_1^{⊺}\mathbf{x})d\mathbf{x}\\ & =log{f}_p(A_p({\kappa }_0){\mathbf{\mu }}_0,{\mathbf{\mu }}_0,{\kappa }_0)-log{C}_p({\kappa }_1)-{\kappa }_1{\mathbf{\mu }}_1^{⊺}{\int }_{球}\mathbf{x}{f}_p(\mathbf{x};{\mathbf{\mu }}_0,{\kappa }_0)d\mathbf{x}\\ & =log{f}_p(A_p({\kappa }_0){\mathbf{\mu }}_0,{\mathbf{\mu }}_0,{\kappa }_0)-log{C}_p({\kappa }_1)-{\kappa }_1{\mathbf{\mu }}_1^{⊺}{A}_p({\kappa }_0){\mathbf{\mu }}_0\\ & =log{f}_p(A_p({\kappa }_0){\mathbf{\mu }}_0,{\mathbf{\mu }}_0,{\kappa }_0)-log[C_p({\kappa }_1)exp({\kappa }_1{\mathbf{\mu }}_1^{⊺}{A}_p({\kappa }_0){\mathbf{\mu }}_0)]\\ & =log{f}_p(A_p({\kappa }_0){\mathbf{\mu }}_0,{\mathbf{\mu }}_0,{\kappa }_0)-log{f}_p(A_p({\kappa }_0){\mathbf{\mu }}_0,{\mathbf{\mu }}_1,{\kappa }_1)\\ & =log\frac{f_p(A_p({\kappa }_0){\mathbf{\mu }}_0,{\mathbf{\mu }}_0,{\kappa }_0)}{f_p(A_p({\kappa }_0){\mathbf{\mu }}_0,{\mathbf{\mu }}_1,{\kappa }_1)}\end{array}$$ 小结: $\{\begin{array}{l}⟨-log{f}_p(\mathbf{x};\mathbf{\mu },\kappa ){⟩}_{\mathbf{x}\sim vMF(\mathbf{\mu },\kappa )}=-log{f}_p(A_p(\kappa )\mathbf{\mu },\mathbf{\mu },\kappa )\\ {⟨log\frac{f_p(\mathbf{x};{\mathbf{\mu }}_0,{\kappa }_0)}{f_p(\mathbf{x};\mathbf{\mu },{\kappa }_1)}⟩}_{\mathbf{x}\sim vMF({\mathbf{\mu }}_0,{\kappa }_0)}=log\frac{f_p(A_p({\kappa }_0){\mathbf{\mu }}_0,{\mathbf{\mu }}_0,{\kappa }_0)}{f_p(A_p({\kappa }_0){\mathbf{\mu }}_0,{\mathbf{\mu }}_1,{\kappa }_1)}\end{array}$

5. 采样

直接对 $vMF$ 分布进行采样几乎是不可能的. 先看 Wikipedia 怎么说:

大概意思是说, 使用一个叫"radial-tangential decomposition"的东西, 它将 $\mathbf{x}$ 分解为 $\mathbf{x}=t\mathbf{\mu }+\sqrt{1-t^2}\mathbf{v}$, 其中 $\mathbf{v}\in {\mathbb{R}}^p$ 位于 $(p-2)$ 维的单位切子球上, 这个单位切子球以 $\mathbf{\mu }$ 为中心, 且垂直于 $\mathbf{\mu }$. 从 $[-1,1]$ 按概率密度 $f_{radial}(t;\kappa ,p)$ 取一个 $t$、切子球上均匀地取一个 $\mathbf{v}$, 就完成了 $\mathbf{x}$ 的采样.

这太过抽象, 让人一头雾水, 好在从 Fast Python Sampler of the von Mises Fisher Distribution 中找到了一张图:

所谓的切子球就是图中绿色圆圈和球相交的圆, 蓝色的向量为 $t\mathbf{\mu }$, 黄色的向量为 $\sqrt{1-t^2}\mathbf{v}$. 可见 $\mathbf{v}$ 确实在切子球上的均匀分布. 若想对 $vMF$ 分布采样的话, 主要在于 $t$ 的采样和相应的切子球上的均匀采样.

根据 Wikipedia, $t$ 的概率密度函数是: $$\begin{array}{rlr}{f}_{radial}(t;\kappa ,p)& =\frac{(\kappa /2{)}^{\nu }}{\Gamma (\frac{1}{2})\Gamma (\nu +\frac{1}{2})I_{\nu }(\kappa )}(1-t^2{)}^{\nu -\frac{1}{2}}exp(\kappa t)& \nu =\frac{p}{2}-1\\ & =\frac{(\kappa /2{)}^{\frac{p}{2}-1}}{\sqrt{\pi }\Gamma (\frac{p-1}{2})I_{\frac{p}{2}-1}(\kappa )}(1-t^2{)}^{\frac{p-3}{2}}exp(\kappa t)& \end{array}$$ 比较复杂, 对于 $t$ 的采样, 下面会单独详细地讲解, 这里暂时先假设已经采样到 $t$ 样本了.

然而, 对于一般的 $\mathbf{\mu }$, 式子 $\mathbf{x}=t\mathbf{\mu }+\sqrt{1-t^2}\mathbf{v}$ 的计算并不容易, 因为 $\mathbf{v}$ 的采样比较困难. 但当 $\mathbf{\mu }={\mathbf{e}}_1=(1,0,\cdots ,0)$ 时:

1. 按概率密度 $f_{radial}(t;\kappa ,p)$ 采样一个 $t\in [-1,1]$;
2. 从 $(p-1)$ 元标准正太分布采样一个 $\mathbf{v}$, 并归一化, 就是切子球上的均匀分布的子向量;
3. $\mathbf{x}=[t|\sqrt{1-t^2}\mathbf{v}]$ 为所求采样, 其中 $|$ 是拼接.(注意 $\mathbf{v}$ 比 $\mathbf{x}$ 少了一个维度)

然后, 我们可以借助 Householder Transform - Wikipedia 将样本 $\mathbf{x}$ 变换到真正想要的 $vMF(\mathbf{\mu },\kappa )$ 分布. 具体地: $$\begin{array}{rl}{P}_v& =\frac{{\mathbf{e}}_1-\mathbf{\mu }}{\Vert {\mathbf{e}}_1-\mathbf{\mu }\Vert }\\ \mathbf{x}& =\mathbf{x}-2P_v^{⊺}\mathbf{x}{P}_v^{⊺}\end{array}$$ 就是真正想要的样本了. 这相当于将 ${\mathbf{e}}_1$ 转向了 $\mathbf{\mu }$, 其对应的 $vMF({\mathbf{e}}_1,\kappa )$ 也变成了 $vMF(\mathbf{\mu },\kappa )$.

Householder Transform, 这种计算是一种镜像翻转

至此, 只剩一个问题没解决: $t$ 的采样. 由于 $f_{radial}(t;\kappa ,p)$ 的复杂性, 对 $t$ 的采样比较麻烦, 先从特殊情况 $p=3$ 讲起, 由于消除了 $(1-t^2{)}^{\frac{p-3}{2}}$, 利用 Inverse Transform Sampling 可实现高效采样. 但当 $p>3$ 时需要用到 Rejection Sampling, 且相当复杂.

5.1 $p=3$ 时的 Inverse Transform Sampling

当 $p=3$ 时, $t$ 的概率密度函数为: $$\begin{array}{rl}{f}_{radial}(t;\kappa ,3)& =\frac{(\kappa /2{)}^{\frac{1}{2}}}{\sqrt{\pi }\Gamma (1)I_{\frac{1}{2}}(\kappa )}(1-t^2{)}^{0}exp(\kappa t)\\ & =\frac{{\kappa }^{\frac{1}{2}}}{\sqrt{2\pi }{I}_{\frac{1}{2}}(\kappa )}exp(\kappa t)\end{array}$$ 此时 $(1-t^2{)}^{\frac{p-3}{2}}$ 项消失, 可以积出其累积分布函数, 然后利用 Inverse Transform Sampling 就可以高效地完成采样.

Inverse Transform Sampling: 任何连续分布的累积分布函数 $F(x)$ 服从 $U(0,1)$ 分布.

计算其累积分布函数: $$\begin{array}{rl}F(t;\kappa ,3)& ={\int }_{-1}^t{f}_{radial}(x;\kappa ,3)dx\\ & =\frac{{\kappa }^{\frac{1}{2}}}{\sqrt{2\pi }{I}_{\frac{1}{2}}(\kappa )}{\int }_{-1}^t{e}^{\kappa x}dx\\ & =\frac{1}{\sqrt{2\pi \kappa }{I}_{\frac{1}{2}}(\kappa )}{e}^{\kappa x}{|}_{-1}^t\\ & =\frac{e^{\kappa t}-e^{-\kappa }}{\sqrt{2\pi \kappa }{I}_{\frac{1}{2}}(\kappa )}\\ F(1;\kappa ,3)& =\frac{e^{\kappa t}-e^{-\kappa }}{\sqrt{2\pi \kappa }{I}_{\frac{1}{2}}(\kappa )}=1\\ ⟹F(t;\kappa ,3)& =\frac{e^{\kappa t}-e^{-\kappa }}{e^{\kappa }-e^{-\kappa }}\end{array}$$ 按照 Inverse Transform Sampling, 令 $u=F(t;\kappa ,3)\sim U(0,1)$, 得: $$t=\frac{1}{\kappa }ln(u{e}^{\kappa }+(1-u)e^{-\kappa })$$ 与 Wikipedia 中的 是等价的, Wikipedia 的写法不带 $e^{\kappa }$, 当 $\kappa$ 比较大的时候, 也许可以避免溢出.

画出其 $t$ 随 $u$ 变化的图像:

可见, 当 $u\sim U(0,1)$ 时, $t\in (-1,1)$, 且 $\kappa$ 越大, $t$ 越倾向于 $1$. 例如, 当 $\kappa =10$ 时, 随机取一个 $u$, $t$ 几乎必然落在 $[0.5,1]$ 内. 当 $\kappa =0$ 时(可通过洛必达法则算极限), $t=2u-1\sim U(-1,1)$.

5.2 $p>3$ 时的 Rejection Sampling

此时累积分布函数不好求, 只能使用 Rejection Sampling 了.

接受-拒绝法的基本想法如下。假设 $p(x)$ 不可以直接抽样。找一个可以直接抽样的分布，称为建议分布 (proposal distribution)。假设 $q(x)$ 是建议分布的概率密度数，并且有 $q(x)$ 的 $c$ 倍一定大于等于 $p(x)$，其中 $c>0$, 如图所示。

按照 $q(x)$ 进行抽样, 假设得到结果是 $x^{\ast }$, 再按照 $\frac{p(x^{\ast })}{cq(x^{\ast })}$ 的比例随机决定是否接受 $x^{\ast }$。直观上，落到 $p(x^{\ast })$ 范围内的就接受(绿色), 落到 $p(x^{\ast })$ 范围外的就拒绝(红色)。接受-拒绝法实际是按照 $p(x)$ 的涵盖面积(或涵盖体积)占 $cq(x)$ 的涵盖面积(或涵盖体积)的比例进行抽样。
假设有 $x_1,x_2$, 其实我们需要的是: 采到 $x_1$ 的概率和采到 $x_2$ 的概率比值是 $\frac{p(x_1)}{p(x_2)}$, 我们来看一看按照上面的方式采样, 到底是不是这样. 首先, 根据 $q(x)$ 采样得到 $x_1$ 的概率是 $q(x_1)$, 被接受的概率是 $\frac{p(x_1)}{cq(x_1)}$, 那么整体来讲采到 $x_1$ 的概率是 $q(x_1)\frac{p(x_1)}{cq(x_1)}=\frac{p(x_1)}{c}$, 同理得到 $x_2$ 的概率是 $\frac{p(x_2)}{c}$, 那么两者比值为 $\frac{p(x_1)}{c}/\frac{p(x_2)}{c}=\frac{p(x_1)}{p(x_2)}$.
接受率的均值: $$\int \frac{p(x)}{cq(x)}q(x)dx=\frac{1}{c}$$ 所以, 在保证 $cq(x)\ge p(x)$ 的同时, 要尽量使 $c$ 小, 以使接受概率最大化.

下面介绍来自论文 Computer Generation of Distributions on the m-spher 的接受-拒绝采样法. 其所选的 proposal distribution 概率密度函数为: $$\begin{array}{rl}e(t;b)=& \frac{2b^{\frac{p-1}{2}}}{\mathrm{B}(\frac{p-1}{2},\frac{p-1}{2})}\frac{(1-t^2{)}^{\frac{p-3}{2}}}{[(1+b)-(1-b)t{]}^{p-1}}其中t\in (-1,1),b\in (0,1)\end{array}$$ 论文说: 从 $Beta$ 分布 $Be(\frac{p-1}{2},\frac{p-1}{2})$, 中采样得到 $y$, 再计算 $t=\frac{1-(1+b)y}{1-(1-b)y}$, 则 $t$ 满足概率密度函数 $e(x,b)$. 至于 $Beta$ 分布的采样, 有不少 Python 包可以使用, 所以, “可以直接抽样的分布” 找到了.

另外, 为了最大化接受率, 求得最佳 $b_0=\frac{-2\kappa +\sqrt{4{\kappa }^2+(p-1{)}^2}}{p-1}$, 经过一番复杂的计算, 论文给出的采样算法是:

1. Initialization: Set $$\begin{array}{rl}b& =\frac{-2\kappa +\sqrt{4{\kappa }^2+(p-1{)}^2}}{p-1}\\ a& =\frac{(p-1)+2\kappa +\sqrt{4{\kappa }^2+(p-1{)}^2}}{4}\\ d& =\frac{4ab}{(1+b)-(p-1)log(p-1)}\end{array}$$
2. Generate $y\sim Be(\frac{p-1}{2},\frac{p-1}{2})$ and $u\sim U(0,1)$. Compute $$t=\frac{1-(1+b)y}{1-(1-b)y}ands=\frac{2ab}{1+(1+b)y}$$
3. If $(p-1)logs-s+d<logu$, go to step 2. else accept $t$.

经过代码试验, 算法是对的, 包括:

• 采 1000000 个样本, 其均值满足 $mean=A_p(\kappa )\mathbf{\mu }$;
• 计算这 1000000 个样本与 $\mathbf{\mu }$ 的内积, 并计算内积的均值, 再使用 scipy.stats.vonmises_fisher 进行同样的操作, 发现内积均值相等, 说明两者分布应该是一样的.

详细推导见 $vMF$ 分布的接受-拒绝采样, 其中还包括拒绝采样的平均接受率.

5.3 Fast Python Sampler of the von Mises Fisher Distribution

对于 $t$ 的采样, 论文 Fast Python Sampler of the von Mises Fisher Distribution 提出了另一种拒绝采样方法:

大致原理是一样的, 也是将分布 $Be(\frac{p-1}{2},\frac{p-1}{2})$ 变换后, 作为建议分布, 执行拒绝采样, 论文声称速度会更快, 但经过推导, 发现其和上面的采样算法是一模一样的. 其详细推导见 Appendix 5.3.