von Mises-Fisher Distribution (Appendix 1)

2. Relation to Normal Distribution

疑问：有没有不各向同性的 vMF？
答：应该是没有的，如果想让各方向偏离中心的速度不一致，则协方差矩阵不为 $\mathbf{I}$ 的倍数. 正态分布的概率密度函数为：$$\begin{array}{r}f(\mathbf{x})=\frac{1}{\sqrt{(2\pi {)}^p|\Sigma |}}{e}^{-\frac{1}{2}(\mathbf{x}-\mathbf{\mu }{)}^{⊺}{\Sigma }^{-1}(\mathbf{x}-\mathbf{\mu })}\end{array}$$ 类比上面的推导，我们需要得出形似：$$\begin{array}{rl}{G}_p(\mathbf{x};\mathbf{\mu },\kappa )& =C(p,\kappa ,r)exp\left(\kappa r\frac{{\mathbf{\mu }}^{⊺}}{r}{\Sigma }^{-1}\mathbf{x}\right)\end{array}$$ 的东西，所以，必要的，需要：$$\begin{array}{rl}{\mathbf{x}}^{⊺}{\Sigma }^{-1}\mathbf{x}& =const\end{array}$$ 我们知道，要想让 ${\mathbf{x}}^{⊺}{\Sigma }^{-1}\mathbf{x}=const$ 表示为球，则必须使 ${\Sigma }^{-1}=\alpha \mathbf{I}$，所以，假设不了“协方差矩阵不为 $\mathbf{I}$ 的倍”, 也就不可能存在不各向同性的 vMF.

[注]：${\mathbf{x}}^{⊺}{\Sigma }^{-1}\mathbf{x}=const$，表示一般的超椭球。

等等！对于 $\mathbf{x}$ 来说, 即使它不是单位向量，也代表了一个方向，$\frac{\mathbf{x}}{\Vert \mathbf{x}\Vert }$ 的分布会是非各向同性的 vMF 吗？有可能是哎！
假设 $\mathbf{y}=\left[\begin{array}{c}{y}_1\\ y_2\end{array}\right]=\left[\begin{array}{cc}1& 0\\ 0& 2\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]=\mathbf{Ax}$，$\mathbf{y}$ 代表单位圆的话，则 $\mathbf{x}$ 在椭圆 $x_1^2+(2x_2{)}^2=1$ 上，如下图：

可以看到，B 和 D 点是对应的，那么弧 AB 和 AD 的“点数”应该是一样的，而弧 AB 对应的方向却在弧 AC 上，即 AD 段对应的方向，被压缩在了 AC 段，如果采样 $\mathbf{y}$ 的话，对应的 $\mathbf{x}$ 方向会更集中趋于椭圆的长轴，实现了“非各向同性”，只不过分布在椭圆上，且采样后需要归一化处理。

5. 采样

$t$ 的概率密度函数推导

直接将 $t={\mathbf{\mu }}^{⊺}\mathbf{x}$ 代入 $f_p(\mathbf{x};\mathbf{\mu },\kappa )$, 得: $$\begin{array}{rlr}{f}_p(\mathbf{x};\mathbf{\mu },\kappa )& =C_p(\kappa )exp(\kappa {\mathbf{\mu }}^{⊺}\mathbf{x})& \\ & =C_p(\kappa )exp(\kappa t)& t\in [-1,1]\\ & =C_p(\kappa )exp(\kappa cos\theta )& \theta \in [0,\pi ]\end{array}$$ 注意这是 $\mathbf{x}$ 一点处的概率密度. 沿着 $t$ 处的切子球求积分, 以得到 $t$ 或 $\theta$ 处的整个概率密度: $$\begin{array}{rlr}{\int }_{切子球}{f}_p(\mathbf{x};\mathbf{\mu },\kappa )ds& ={\int }_{切子球}{C}_p(\kappa )exp(\kappa \mathbf{x})ds& \\ & =C_p(\kappa )exp(\kappa t)\frac{2{\pi }^{\frac{p-1}{2}}}{\Gamma (\frac{p-1}{2})}(1-t^2{)}^{\frac{p-2}{2}}& S^{p-2}的表面积\propto r^{p-2}\\ & =C_p(\kappa )exp(\kappa cos\theta )\frac{2{\pi }^{\frac{p-1}{2}}}{\Gamma (\frac{p-1}{2})}si{n}^{p-2}\theta & \end{array}$$ 根据 n-sphere - Wikipedia, 切子球 $S^{p-2}$ 的表面积 $S_{p-2}=\frac{2{\pi }^{\frac{p-1}{2}}}{\Gamma (\frac{p-1}{2})}{r}^{p-2}$, 再沿 $t$ 或 $\theta$ 积分: $$\begin{array}{rlr}& {\int }_0^{\pi }{C}_p(\kappa )exp(\kappa cos\theta )\frac{2{\pi }^{\frac{p-1}{2}}}{\Gamma (\frac{p-1}{2})}si{n}^{p-2}\theta d\theta & \\ =& C_p(\kappa )\frac{2{\pi }^{\frac{p-1}{2}}}{\Gamma (\frac{p-1}{2})}{\int }_1^{-1}exp(\kappa t)(1-t^2{)}^{\frac{p-2}{2}}(\frac{-1}{\sqrt{1-t^2}}dt)& \because cos0=1,cos\pi =-1\\ =& C_p(\kappa )\frac{2{\pi }^{\frac{p-1}{2}}}{\Gamma (\frac{p-1}{2})}{\int }_{-1}^{1}exp(\kappa t)(1-t^2{)}^{\frac{p-3}{2}}dt& \\ =& C_p(\kappa )\frac{2{\pi }^{\nu +\frac{1}{2}}}{\Gamma (\nu +\frac{1}{2})}{\int }_{-1}^{1}exp(\kappa t)(1-t^2{)}^{\nu -\frac{1}{2}}dt& 令\nu =\frac{p}{2}-1\end{array}$$ 那么, 将 $$\begin{array}{rlr}{\int }_{-1}^{1}exp(\kappa t)(1-t^2{)}^{\nu -\frac{1}{2}}dt& =\frac{\Gamma (\frac{1}{2})\Gamma (\nu +\frac{1}{2})I_{\nu }(\kappa )}{(\kappa /2{)}^{\nu }}& [I_{\nu }的公式]\\ & =\frac{\Gamma (\frac{1}{2})\Gamma (\frac{p-1}{2})I_{\frac{p}{2}-1}(\kappa )}{(\kappa /2{)}^{\frac{p}{2}-1}}& \end{array}$$ 代入, 得: $$\begin{array}{rlr}& C_p(\kappa )\frac{2{\pi }^{\nu +\frac{1}{2}}}{\Gamma (\nu +\frac{1}{2})}{\int }_{-1}^{1}exp(\kappa t)(1-t^2{)}^{\nu -\frac{1}{2}}dt& \\ =& C_p(\kappa )\frac{2{\pi }^{\nu +\frac{1}{2}}}{\Gamma (\nu +\frac{1}{2})}\frac{\Gamma (\frac{1}{2})\Gamma (\nu +\frac{1}{2})I_{\nu }(\kappa )}{(\kappa /2{)}^{\nu }}& \\ =& \frac{2{\pi }^{\nu +1}{I}_{\nu }(\kappa )}{(\kappa /2{)}^{\nu }}{C}_p(\kappa )& \Gamma (\frac{1}{2})={\pi }^{\frac{1}{2}}\\ =& \frac{2{\pi }^{\frac{p}{2}}{I}_{\frac{p}{2}-1}(\kappa )}{(\kappa /2{)}^{\frac{p}{2}-1}}{C}_p(\kappa )& \\ =& \frac{(2\pi {)}^{\frac{p}{2}}{I}_{\frac{p}{2}-1}(\kappa )}{{\kappa }^{\frac{p}{2}-1}}{C}_p(\kappa )=1& (2)\\ ⟺& & \\ C_p(\kappa )=& \frac{{\kappa }^{\frac{p}{2}-1}}{(2\pi {)}^{\frac{p}{2}}{I}_{\frac{p}{2}-1}(\kappa )}& [与Wikipedia一致]\end{array}$$ 所以, 由式 $(2)$ 可得 $t$ 的概率密度函数为: $$\begin{array}{rlr}{f}_{radial}(t;\kappa ,p)& =C_p(\kappa )\frac{2{\pi }^{\nu +\frac{1}{2}}}{\Gamma (\nu +\frac{1}{2})}exp(\kappa t)(1-t^2{)}^{\nu -\frac{1}{2}}& \nu =\frac{p}{2}-1\\ & =\frac{(\kappa /2{)}^{\nu }}{\Gamma (\frac{1}{2})\Gamma (\nu +\frac{1}{2})I_{\nu }(\kappa )}exp(\kappa t)(1-t^2{)}^{\nu -\frac{1}{2}}& [与Wikipedia一致]\\ & =\frac{(\kappa /2{)}^{\frac{p}{2}-1}}{\sqrt{\pi }\Gamma (\frac{p-1}{2})I_{\frac{p}{2}-1}(\kappa )}exp(\kappa t)(1-t^2{)}^{\frac{p-3}{2}}& \end{array}$$

单位球上的均匀分布采样

实际上, 从标准正态分布 $Normal(0,1)$ 中采样 $p$ 个数, 组成向量 $\mathbf{v}$, 归一化, 就实现了球 $S_{p-1}$ 上的均匀采样, 并不需要多元正态分布 $MultivariateNormal(0,{\mathbf{I}}^{(p)})$. 因为 $p$ 个独立采样自 $Normal(0,1)$ 的 $p$ 个数概率是: $$\begin{array}{rl}& \prod _{i=1}^p\frac{1}{\sqrt{2\pi }}exp(-\frac{x_i^2}{2})\\ =& \frac{1}{(2\pi {)}^{\frac{p}{2}}}exp\left(-\frac{{\sum }_{i=1}^p{x}_i^2}{2}\right)\end{array}$$ 恰好是 $p$ 元标准正态分布的概率密度.

5.2 拒绝采样

为了与原论文一致, 本小节使用 $m$ 代替 $p$.

首先, 设 $Y\sim Be(\alpha ,\beta )$, 概率密度函数为: $$\begin{array}{rl}f(y;\alpha ,\beta )& =\frac{y^{\alpha -1}(1-y{)}^{\beta -1}}{{\int }_0^1{u}^{\alpha -1}(1-u{)}^{\beta -1}du}\\ & =\frac{\Gamma (\alpha +\beta )}{\Gamma (\alpha )\Gamma (\beta )}{y}^{\alpha -1}(1-y{)}^{\beta -1}\\ & =\frac{1}{\mathrm{B}(\alpha ,\beta )}{y}^{\alpha -1}(1-y{)}^{\beta -1}\end{array}$$ 其中 $y\in (0,1),\alpha ,\beta >0$. 令 $$\begin{array}{rl}X=\frac{1-(1+b)Y}{1-(1-b)Y}& \iff Y=\frac{1-X}{(1+b)-(1-b)X}\\ b& \in (0,1)\end{array}$$ 代入上述概率密度函数 $f(y;\alpha ,\beta )$, 并令 $\alpha =\beta =\frac{m-1}{2}$ 得: $$\begin{array}{rl}& f(\frac{1-x}{(1+b)-(1-b)x};\frac{m-1}{2},\frac{m-1}{2})\\ =& \frac{1}{\mathrm{B}(\frac{m-1}{2},\frac{m-1}{2})}{\left(\left(\frac{1-x}{(1+b)-(1-b)x}\right)\ast \left(1-\frac{1-x}{(1+b)-(1-b)x}\right)\right)}^{\frac{m-1}{2}-1}\\ =& \frac{1}{\mathrm{B}(\frac{m-1}{2},\frac{m-1}{2})}{\left(\left(\frac{1-x}{(1+b)-(1-b)x}\right)\ast \left(\frac{(1+b)-(1-b)x-(1-x)}{(1+b)-(1-b)x}\right)\right)}^{\frac{m-3}{2}}\\ =& \frac{1}{\mathrm{B}(\frac{m-1}{2},\frac{m-1}{2})}{\left(\left(\frac{1-x}{(1+b)-(1-b)x}\right)\ast \left(\frac{b(1+x)}{(1+b)-(1-b)x}\right)\right)}^{\frac{m-3}{2}}\\ =& \frac{1}{\mathrm{B}(\frac{m-1}{2},\frac{m-1}{2})}{\left(\frac{b(1-x^2)}{[(1+b)-(1-b)x{]}^2}\right)}^{\frac{m-3}{2}}\\ =& \frac{b^{\frac{m-3}{2}}}{\mathrm{B}(\frac{m-1}{2},\frac{m-1}{2})}\frac{(1-x^2{)}^{\frac{m-3}{2}}}{[(1+b)-(1-b)x{]}^{m-3}}\\ =& e(x,b)其中x\in (-1,1),b\in (0,1)\end{array}$$ 与论文 Computer Generation of Distributions on the m-spher [1] 不一致, 但我不知道问题出在哪里. 论文里明明说:

却给出了:

由拒绝采样法的 $f_{radial}(x;\kappa ,m)\le M\ast e(x,b)$, 计算: $$\begin{array}{rl}M& =\underset{x}{\max }\frac{f_{radial}(x;\kappa ,m)}{e(x,b)}\\ & =\underset{x}{\max }\frac{\frac{(\kappa /2{)}^{\nu }}{\Gamma (\frac{1}{2})\Gamma (\nu +\frac{1}{2})I_{\nu }(\kappa )}(1-x^2{)}^{\frac{m-3}{2}}exp(\kappa x)}{\frac{b^{\frac{m-3}{2}}}{\mathrm{B}(\frac{m-1}{2},\frac{m-1}{2})}\frac{(1-x^2{)}^{\frac{m-3}{2}}}{[(1+b)-(1-b)x{]}^{m-3}}}\\ & =\underset{x}{\max }\frac{\frac{(\kappa /2{)}^{\nu }}{\Gamma (\frac{1}{2})\Gamma (\nu +\frac{1}{2})I_{\nu }(\kappa )}exp(\kappa x)}{\frac{b^{\frac{m-3}{2}}}{\mathrm{B}(\frac{m-1}{2},\frac{m-1}{2})}\frac{1}{[(1+b)-(1-b)x{]}^{m-3}}}\\ & =\underset{x}{\max }\frac{(\kappa /2{)}^{\nu }}{\Gamma (\frac{1}{2})\Gamma (\nu +\frac{1}{2})I_{\nu }(\kappa )}\frac{\mathrm{B}(\frac{m-1}{2},\frac{m-1}{2})}{b^{\frac{m-3}{2}}}[(1+b)-(1-b)x{]}^{m-3}exp(\kappa x)\end{array}$$ 接下来求极值点: $$\begin{array}{rl}& \frac{\partial \left([(1+b)-(1-b)x{]}^{m-3}exp(\kappa x)\right)}{\partial x}\\ =& (m-3)[(1+b)-(1-b)x{]}^{m-4}\ast [-(1-b)]exp(\kappa x)\\ +& [(1+b)-(1-b)x{]}^{m-3}\kappa exp(\kappa x)=0\\ \Rightarrow & \\ & (m-3)\ast [-(1-b)]+[(1+b)-(1-b)x]\kappa =0\\ & (1+b)-(1-b)x=\frac{(m-3)(1-b)}{\kappa }\\ & (1-b)x=(1+b)-\frac{(m-3)(1-b)}{\kappa }\\ & 解:x_0=\frac{1+b}{1-b}-\frac{m-3}{\kappa }\end{array}$$ 因为 $b\in (0,1),x\in [-1,1]$, 所以 $(1+b)-(1-b)x\ge 0$, 可知, 在经过简约之后, $\frac{\partial \left([(1+b)-(1-b)x{]}^{m-3}exp(\kappa x)\right)}{\partial x}$ 的符号不变, $(m-3)\ast [-(1-b)]+[(1+b)-(1-b)x]\kappa$ 是线性的, 单调递减. 那么 $\frac{f_{radial}(x;\kappa ,m)}{e(x,b)}$ 先曾后减, 在 $x_0$ 处取得最大值.

代入之后, 求得 $M$, 为使接受率最大化(最小化$M$), 对 $b$ 求导: $$\begin{array}{rl}& \frac{\partial \left(\frac{[(1+b)-(1-b)x_0{]}^{m-3}exp(\kappa x_0)}{b^{\frac{m-3}{2}}}\right)}{\partial b}\\ =& exp(-\frac{m-3}{\kappa })\frac{\partial \left(\frac{[\frac{(m-3)(1-b)}{\kappa }{]}^{m-3}exp(\kappa \frac{1+b}{1-b})}{b^{\frac{m-3}{2}}}\right)}{\partial b}\\ =& {\left(\frac{m-3}{\kappa }\right)}^{m-3}exp(-\frac{m-3}{\kappa })\frac{\partial \left(\frac{(1-b{)}^{m-3}exp(\kappa \frac{1+b}{1-b})}{b^{\frac{m-3}{2}}}\right)}{\partial b}=0\\ \Rightarrow & \\ =& \frac{\partial \left(\frac{(1-b{)}^{m-3}exp(\kappa \frac{1+b}{1-b})}{b^{\frac{m-3}{2}}}\right)}{\partial b}=0\\ =& \frac{m-3}{2}(1-\frac{1}{b^2})(b+\frac{1}{b}-2{)}^{\frac{m-5}{2}}exp(\kappa \frac{1+b}{1-b})+(b+\frac{1}{b}-2{)}^{\frac{m-3}{2}}\frac{2\kappa }{(1-b{)}^2}exp(\kappa \frac{1+b}{1-b})\\ \Rightarrow & \\ & \frac{m-3}{2}(1-\frac{1}{b^2})+(b+\frac{1}{b}-2)\frac{2\kappa }{(1-b{)}^2}=0\\ & (m-3)(1-\frac{1}{b^2})+\frac{4\kappa }{b}=0\\ & (m-3)(b^2-1)+4\kappa b=0\\ =& (m-3)b^2+4\kappa b-(m-3)\\ & 解:b_0=\frac{-2\kappa +\sqrt{4{\kappa }^2+(m-3{)}^2}}{m-3}\in (0,1)\end{array}$$

如果把 $(m-3)>0$ 和 $2\kappa >0$ 都看作三角形的直角边, $\sqrt{4{\kappa }^2+(m-3{)}^2}$ 为斜边, 则 $b_0=sec\theta -tan\theta$, 又 $\theta \in (0,\frac{\pi }{2})$, 故 $b_0\in (0,1)$

所有的简约都不改变 $\frac{\partial M}{\partial b}$ 的符号, 所以, $M$ 随 b 先减后增, $b_0$ 是 $M$ 的最小值点.

以上所有的计算中, 对于原论文中的 $e(x,b)$ 和本博文计算的 $e(x,b)$, 差别只在于 $(m-1)$ 和 $(m-3)$, 其他都一样. 继续计算 $x_0$: $$\begin{array}{rl}1+b_0& =1+\frac{-2\kappa +\sqrt{4{\kappa }^2+(m-3{)}^2}}{m-3}\\ & =\frac{m-3-2\kappa +\sqrt{4{\kappa }^2+(m-3{)}^2}}{m-3}\\ & =\frac{m-3+\sqrt{4{\kappa }^2+(m-3{)}^2}-2\kappa }{m-3}\\ 1-b_0& =1-\frac{-2\kappa +\sqrt{4{\kappa }^2+(m-3{)}^2}}{m-3}\\ & =\frac{m-3+2\kappa -\sqrt{4{\kappa }^2+(m-3{)}^2}}{m-3}\\ & =\frac{m-3+2\kappa -\sqrt{4{\kappa }^2+(m-3{)}^2}}{m-3}\\ x_0& =\frac{1+b_0}{1-b_0}-\frac{m-3}{\kappa }\\ & =\frac{m-3-2\kappa +\sqrt{4{\kappa }^2+(m-3{)}^2}}{m-3+2\kappa -\sqrt{4{\kappa }^2+(m-3{)}^2}}-\frac{m-3}{\kappa }\\ & =\frac{(m-3+\sqrt{4{\kappa }^2+(m-3{)}^2}-2\kappa )(m-3+\sqrt{4{\kappa }^2+(m-3{)}^2}+2\kappa )}{(m-3+2\kappa {)}^2-(4{\kappa }^2+(m-3{)}^2)}-\frac{m-3}{\kappa }\\ & =\frac{(m-3+\sqrt{4{\kappa }^2+(m-3{)}^2}{)}^2-4{\kappa }^2}{(m-3{)}^2+4{\kappa }^2+4\kappa (m-3)-(4{\kappa }^2+(m-3{)}^2)}-\frac{m-3}{\kappa }\\ & =\frac{(m-3{)}^2+4{\kappa }^2+(m-3{)}^2+2(m-3)\sqrt{4{\kappa }^2+(m-3{)}^2}-4{\kappa }^2}{4\kappa (m-3)}-\frac{m-3}{\kappa }\\ & =\frac{2(m-3{)}^2+2(m-3)\sqrt{4{\kappa }^2+(m-3{)}^2}}{4\kappa (m-3)}-\frac{m-3}{\kappa }\\ & =\frac{(m-3)+\sqrt{4{\kappa }^2+(m-3{)}^2}}{2\kappa }-\frac{2(m-3)}{2\kappa }\\ & =\frac{-(m-3)+\sqrt{4{\kappa }^2+(m-3{)}^2}}{2\kappa }\in (0,1)\\ \frac{1-b_0}{1+b_0}=& \frac{m-3+2\kappa -\sqrt{4{\kappa }^2+(m-3{)}^2}}{m-3-2\kappa +\sqrt{4{\kappa }^2+(m-3{)}^2}}\\ =& \frac{[m-3+2\kappa -\sqrt{4{\kappa }^2+(m-3{)}^2}][(m-3-2\kappa )-\sqrt{4{\kappa }^2+(m-3{)}^2}]}{[(m-3-2\kappa )+\sqrt{4{\kappa }^2+(m-3{)}^2}][(m-3-2\kappa )-\sqrt{4{\kappa }^2+(m-3{)}^2}]}\\ =& \frac{[(m-3-\sqrt{4{\kappa }^2+(m-3{)}^2})+2\kappa ][(m-3-\sqrt{4{\kappa }^2+(m-3{)}^2})-2\kappa ]}{(m-3-2\kappa {)}^2-[4{\kappa }^2+(m-3{)}^2]}\\ =& \frac{(m-3-\sqrt{4{\kappa }^2+(m-3{)}^2}{)}^2-4{\kappa }^2}{[(m-3{)}^2+4{\kappa }^2]-4\kappa (m-3)-[4{\kappa }^2+(m-3{)}^2]}\\ =& \frac{(m-3{)}^2+4{\kappa }^2+(m-3{)}^2-2(m-3)\sqrt{4{\kappa }^2+(m-3{)}^2}-4{\kappa }^2}{-4\kappa (m-3)}\\ =& \frac{2(m-3{)}^2-2(m-3)\sqrt{4{\kappa }^2+(m-3{)}^2}}{-4\kappa (m-3)}\\ =& \frac{-(m-3)+\sqrt{4{\kappa }^2+(m-3{)}^2}}{2\kappa }=x_0\end{array}$$ 采样一个 $u\sim Uniform(0,1)$, $y\sim Be(\frac{m-1}{2},\frac{m-1}{2})\Rightarrow x=\frac{1-(1+b_0)y}{1-(1-b_0)y}\sim e(x,b_0)$ 当 $u<\frac{f_{radial}(x;\kappa ,m)}{M\ast e(x,b_0)}$ 时, 接受 $x$ 为样本. $$\begin{array}{rl}\frac{f_{radial}(x;\kappa ,m)}{M\ast e(x,b_0)}& =\frac{f_{radial}(x;\kappa ,m)}{Mf(\frac{1-x}{(1+b_0)-(1-b_0)x};\frac{m-1}{2},\frac{m-1}{2})}\\ & =\frac{(\kappa /2{)}^{\nu }}{\Gamma (\frac{1}{2})\Gamma (\nu +\frac{1}{2})I_{\nu }(\kappa )}\frac{\mathrm{B}(\frac{m-1}{2},\frac{m-1}{2})}{b_0^{\frac{m-3}{2}}}[(1+b_0)-(1-b_0)x{]}^{m-3}exp(\kappa x)/M\\ & =\frac{[(1+b_0)-(1-b_0)x{]}^{m-3}exp(\kappa x)}{[(1+b_0)-(1-b_0)x_0{]}^{m-3}exp(\kappa x_0)}\\ & ={\left(\frac{(1+b_0)-(1-b_0)x}{(1+b_0)-(1-b_0)x_0}\right)}^{m-3}exp(\kappa (x-x_0))\\ log\frac{f_{radial}(x;\kappa ,m)}{M\ast e(x,b_0)}& =(m-3)log\left(\frac{(1+b_0)-(1-b_0)x}{(1+b_0)-(1-b_0)x_0}\right)+\kappa (x-x_0)\\ & =(m-3)log\left(\frac{1-\frac{1-b_0}{1+b_0}x}{1-\frac{1-b_0}{1+b_0}{x}_0}\right)+\kappa (x-x_0)\\ & =(m-3)log\left(\frac{1-x_{0}x}{1-x_0^2}\right)+\kappa (x-x_0)\\ & =(m-3)log(1-x_{0}x)-(m-3)log(1-x_0^2)+\kappa (x-x_0)\\ & =\kappa x+(m-3)log(1-x_{0}x)-[\kappa x_0+(m-3)log(1-x_0^2)]\\ & =\kappa x+(m-3)log(1-x_{0}x)-c\end{array}$$ 即, 当 $logu<\kappa x+(m-3)log(1-x_{0}x)-c$ 时接受样本 $x$. 与 Communications in Statistics - Simulation and Computation [2] 中的判别公式一致了(除了 $(m-3)$ 和 $(m-1)$):

但是与 [1] 不一致, 且 [2] 声称 [1] 的采样方案有问题:
意思大概是作者在实践中并没有的到正确的采样, 且不知道问题出在哪里, 所以他们又根据 [1] 推导了一遍, 得出了上述接受-拒绝判别式. 但, 我也推导了一遍, 发现两者是等价的…

问题到底在哪里?

似乎所有的差别就在 $e(x,b)$, 到底是: $$\begin{array}{rl}e(x,b)=& \frac{b^{\frac{m-3}{2}}}{\mathrm{B}(\frac{m-1}{2},\frac{m-1}{2})}\frac{(1-x^2{)}^{\frac{m-3}{2}}}{[(1+b)-(1-b)x{]}^{m-3}}\end{array}$$ 还是: $$\begin{array}{rl}e(x,b)=& \frac{2b^{\frac{m-1}{2}}}{\mathrm{B}(\frac{m-1}{2},\frac{m-1}{2})}\frac{(1-x^2{)}^{\frac{m-3}{2}}}{[(1+b)-(1-b)x{]}^{m-1}}\end{array}$$ ? ? ? ? ? ? ? 找到问题出在哪了, 保留以上过程, 以警钟长鸣 ? ? ? ? ? ? ? ?
记住, 永远不要想当然! 概率密度函数的变换不是进行一下变量代换就 OK 了的! 前面 $t$ 的概率密度函数推导中, 就知道 $t$ 和 $\theta$ 在变换的时候需要带上 $dt,d\theta$, 怎么现在就忘了? 根本原因是对知识的理解不够透彻. 我们的研究对象本质上是概率 $f(y;\alpha ,\beta )dy$, 而不是概率密度 $f(y;\alpha ,\beta )$, 当进行变量代换时, 不光要变换概率密度中的变量, $dy$ 也要变换, 甚至求 $[a,b]$ 区间上的概率时, ${\int }_a^{b}f(y;\alpha ,\beta )dy$, 积分区间也要变换.

在得到随机变量变换关系: $$\begin{array}{rl}X=\frac{1-(1+b)Y}{1-(1-b)Y}& \iff Y=\frac{1-X}{(1+b)-(1-b)X}\end{array}$$ 后, 假如让你积分求区间 $[y_1,y_2]$ 上的概率呢? $$\begin{array}{rl}{\int }_{y_1}^{y_2}f(y;\alpha ,\beta )dy& ={\int }_{\frac{1-x_1}{(1+b)-(1-b)x_1}}^{\frac{1-x_2}{(1+b)-(1-b)x_2}}f\left(\frac{1-x}{(1+b)-(1-b)x};\alpha ,\beta \right)d\left(\frac{1-x}{(1+b)-(1-b)x}\right)\\ & ={\int }_{\frac{1-x_1}{(1+b)-(1-b)x_1}}^{\frac{1-x_2}{(1+b)-(1-b)x_2}}f\left(\frac{1-x}{(1+b)-(1-b)x};\alpha ,\beta \right)\frac{-2b}{[(1+b)-(1-b)x{]}^2}dx\end{array}$$

$$\begin{array}{rl}dy& =d\left(\frac{1-x}{(1+b)-(1-b)x}\right)\\ & =\frac{-[(1+b)-(1-b)x]-(1-x)[-(1-b)]}{[(1+b)-(1-b)x{]}^2}dx\\ & =\frac{-(1+b)+(1-b)x+(1-b)-(1-b)x}{[(1+b)-(1-b)x{]}^2}dx\\ & =\frac{-2b}{[(1+b)-(1-b)x{]}^2}dx\end{array}$$

由于 $\frac{-2b}{[(1+b)-(1-b)x{]}^2}<0$, 变量 $x,y$ 之间是负相关的, 那么积分式写为: $$\begin{array}{r}{\int }_{\frac{1-x_2}{(1+b)-(1-b)x_2}}^{\frac{1-x_1}{(1+b)-(1-b)x_1}}f\left(\frac{1-x}{(1+b)-(1-b)x};\alpha ,\beta \right)\frac{2b}{[(1+b)-(1-b)x{]}^2}dx\end{array}$$ 所以, 关于 $x$ 的概率密度函数是: $$\begin{array}{rl}e(x,b)=& f\left(\frac{1-x}{(1+b)-(1-b)x};\alpha ,\beta \right)\frac{2b}{[(1+b)-(1-b)x{]}^2}\\ =& \frac{b^{\frac{m-3}{2}}}{\mathrm{B}(\frac{m-1}{2},\frac{m-1}{2})}\frac{(1-x^2{)}^{\frac{m-3}{2}}}{[(1+b)-(1-b)x{]}^{m-3}}\frac{2b}{[(1+b)-(1-b)x{]}^2}\\ =& \frac{2b^{\frac{m-1}{2}}}{\mathrm{B}(\frac{m-1}{2},\frac{m-1}{2})}\frac{(1-x^2{)}^{\frac{m-3}{2}}}{[(1+b)-(1-b)x{]}^{m-1}}\end{array}$$ 这下一致了.

平均接受率

在拒绝采样法的简介中说: 平均接受率为 $\frac{1}{c}$, 那么本拒绝采样法的平均接受率为 $\frac{1}{M}$. $$\begin{array}{r}\frac{1}{M}=\frac{e(x_0,b_0)}{f(x_0,\kappa ,m)}\end{array}$$ 公式太过复杂, 简化无果后, 在 desmos 网站画出了其随 $\kappa ,m$ 变化的图像(计算比较复杂,可能需要等一会儿图像才会出来):

可见, 接受率还是比较高的, 趋势大概是: $m$ 越大, 接受率越高, $\kappa$ 越大, 接受率越低. 但即使在 $m=4,\kappa =20$ 时, 接受率依然大于 $0.7$.

下面给出了 $(m=50,\kappa =5)$ 和 $(m=5,\kappa =20)$ 时, $f_{radial}(x;\kappa ,m)$ 和 $M\ast e(x,b_0)$ 的图像(前者为蓝色, 后者为黑色). 可见, 当 $(m=50,\kappa =5)$ 时, 两个图像几乎是重合的; 当 $(m=5,\kappa =20)$ 时, 两者出现了较大出入, 表现为接受率的下降.

参考文献

[1] Computer Generation of Distributions on the m-spher.
[2] Communications in Statistics - Simulation and Computation.