5.3 Fast Python Sampler of the von Mises Fisher Distribution [3]
从论文中
p
r
o
c
e
d
u
r
e
A
n
g
l
e
G
e
n
e
r
a
t
o
r
(
d
,
κ
)
procedure~AngleGenerator(d, κ)
procedure AngleGenerator(d,κ) 中的变换来看, 假设
y
∼
B
e
(
m
−
1
2
,
m
−
1
2
)
y \sim Be(\frac{m-1}{2}, \frac{m-1}{2})
y∼Be(2m−1,2m−1), 有:
x
=
r
+
(
2
y
−
1
)
1
+
r
(
2
y
−
1
)
⇔
y
=
(
1
−
r
)
(
1
+
x
)
2
(
1
−
r
x
)
\begin{aligned} x = \frac{r+(2y-1)}{1+r(2y-1)} ~\Leftrightarrow~ y = \frac{(1-r)(1+x)}{2(1-rx)} \end{aligned}
x=1+r(2y−1)r+(2y−1) ⇔ y=2(1−rx)(1−r)(1+x) 则
1
−
y
=
1
−
(
1
−
r
)
(
1
+
x
)
2
(
1
−
r
x
)
=
(
1
+
r
)
(
1
−
x
)
2
(
1
−
r
x
)
f
(
y
;
m
−
1
2
,
m
−
1
2
)
d
y
=
1
B
(
m
−
1
2
,
m
−
1
2
)
(
y
(
1
−
y
)
)
m
−
3
2
d
y
=
1
B
(
m
−
1
2
,
m
−
1
2
)
(
(
1
−
r
)
(
1
+
x
)
2
(
1
−
r
x
)
(
1
+
r
)
(
1
−
x
)
2
(
1
−
r
x
)
)
m
−
3
2
d
(
(
1
−
r
)
(
1
+
x
)
2
(
1
−
r
x
)
)
=
1
B
(
m
−
1
2
,
m
−
1
2
)
(
1
−
r
2
)
m
−
3
2
(
1
−
x
2
)
m
−
3
2
[
2
(
1
−
r
x
)
]
m
−
3
2
(
1
−
r
2
)
[
2
(
1
−
r
x
)
]
2
d
x
=
(
1
−
r
2
)
m
−
1
2
2
m
−
2
B
(
m
−
1
2
,
m
−
1
2
)
(
1
−
x
2
)
m
−
3
2
(
1
−
r
x
)
m
−
1
d
x
\begin{aligned} 1-y &= 1 - \frac{(1-r)(1+x)}{2(1-rx)} \\ &= \frac{(1+r)(1-x)}{2(1-rx)} \\ f(y; \frac{m-1}{2}, \frac{m-1}{2})dy &= \frac{1}{B(\frac{m-1}{2}, \frac{m-1}{2})} \left( y(1-y) \right)^{\frac{m-3}{2}} dy \\ &= \frac{1}{B(\frac{m-1}{2}, \frac{m-1}{2})} \left( \frac{(1-r)(1+x)}{2(1-rx)} \frac{(1+r)(1-x)}{2(1-rx)} \right)^{\frac{m-3}{2}} d\left( \frac{(1-r)(1+x)}{2(1-rx)} \right) \\ &= \frac{1}{B(\frac{m-1}{2}, \frac{m-1}{2})} \frac{(1-r^2)^\frac{m-3}{2}(1-x^2)^\frac{m-3}{2}}{[2(1-rx)]^{m-3}} \frac{2(1-r^2)}{[2(1-rx)]^2} dx \\ &= \frac{(1-r^2)^\frac{m-1}{2}}{2^{m-2}B(\frac{m-1}{2}, \frac{m-1}{2})} \frac{(1-x^2)^\frac{m-3}{2}}{(1-rx)^{m-1}} dx \end{aligned}
1−yf(y;2m−1,2m−1)dy=1−2(1−rx)(1−r)(1+x)=2(1−rx)(1+r)(1−x)=B(2m−1,2m−1)1(y(1−y))2m−3dy=B(2m−1,2m−1)1(2(1−rx)(1−r)(1+x)2(1−rx)(1+r)(1−x))2m−3d(2(1−rx)(1−r)(1+x))=B(2m−1,2m−1)1[2(1−rx)]m−3(1−r2)2m−3(1−x2)2m−3[2(1−rx)]22(1−r2)dx=2m−2B(2m−1,2m−1)(1−r2)2m−1(1−rx)m−1(1−x2)2m−3dx 所以, 这里的
e
(
x
,
r
)
e(x, r)
e(x,r) 为:
e
(
x
;
r
)
=
(
1
−
r
2
)
m
−
1
2
2
m
−
2
B
(
m
−
1
2
,
m
−
1
2
)
(
1
−
x
2
)
m
−
3
2
(
1
−
r
x
)
m
−
1
\begin{aligned} e(x; r)= \frac{(1-r^2)^\frac{m-1}{2}}{2^{m-2}B(\frac{m-1}{2}, \frac{m-1}{2})} \frac{(1-x^2)^\frac{m-3}{2}}{(1-rx)^{m-1}} \end{aligned}
e(x;r)=2m−2B(2m−1,2m−1)(1−r2)2m−1(1−rx)m−1(1−x2)2m−3 我们再进行同样的流程:
计算:
M
=
max
x
f
r
a
d
i
a
l
(
x
;
κ
,
m
)
e
(
x
,
b
)
=
max
x
(
κ
/
2
)
ν
Γ
(
1
2
)
Γ
(
ν
+
1
2
)
I
ν
(
κ
)
(
1
−
x
2
)
m
−
3
2
e
x
p
(
κ
x
)
(
1
−
r
2
)
m
−
1
2
2
m
−
2
B
(
m
−
1
2
,
m
−
1
2
)
(
1
−
x
2
)
m
−
3
2
(
1
−
r
x
)
m
−
1
=
max
x
(
κ
/
2
)
ν
Γ
(
1
2
)
Γ
(
ν
+
1
2
)
I
ν
(
κ
)
2
m
−
2
B
(
m
−
1
2
,
m
−
1
2
)
(
1
−
r
2
)
m
−
1
2
(
1
−
r
x
)
m
−
1
e
x
p
(
κ
x
)
\begin{aligned} M &= \max_x \frac{f_{radial}(x;\kappa,m)}{e(x, b)} \\ &= \max_x \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{(1-r^2)^\frac{m-1}{2}}{2^{m-2}B(\frac{m-1}{2}, \frac{m-1}{2})} \frac{(1-x^2)^\frac{m-3}{2}}{(1-rx)^{m-1}} } \\ &= \max_x \frac{(\kappa/2)^\nu}{\Gamma({\frac{1}{2}})\Gamma(\nu+{\frac{1}{2}})I_{\nu}(\kappa)} \frac{2^{m-2}\Beta(\frac{m-1}{2},\frac{m-1}{2})}{(1-r^2)^\frac{m-1}{2}} (1-rx)^{m-1} exp(\kappa x) \end{aligned}
M=xmaxe(x,b)fradial(x;κ,m)=xmax2m−2B(2m−1,2m−1)(1−r2)2m−1(1−rx)m−1(1−x2)2m−3Γ(21)Γ(ν+21)Iν(κ)(κ/2)ν(1−x2)2m−3exp(κx)=xmaxΓ(21)Γ(ν+21)Iν(κ)(κ/2)ν(1−r2)2m−12m−2B(2m−1,2m−1)(1−rx)m−1exp(κx) 接下来求极值点:
∂
[
(
1
−
r
x
)
m
−
1
e
x
p
(
κ
x
)
]
∂
x
=
−
r
(
m
−
1
)
(
1
−
r
x
)
m
−
2
e
x
p
(
κ
x
)
+
(
1
−
r
x
)
m
−
1
κ
e
x
p
(
κ
x
)
=
0
⇒
−
r
(
m
−
1
)
+
(
1
−
r
x
)
κ
=
0
1
−
r
x
=
r
(
m
−
1
)
κ
r
x
=
1
−
r
(
m
−
1
)
κ
解
:
x
0
=
1
r
−
m
−
1
κ
\begin{aligned} & \frac{\partial [(1-rx)^{m-1} exp(\kappa x)]}{\partial x} \\ =& -r(m-1)(1-rx)^{m-2} exp(\kappa x) + (1-rx)^{m-1} \kappa exp(\kappa x) = 0 \\ \Rightarrow \\ & -r(m-1) + (1-rx) \kappa = 0 \\ & 1-rx = \frac{r(m-1)}{\kappa} \\ & rx = 1 - \frac{r(m-1)}{\kappa} \\ & 解:~ x_0 = \frac{1}{r} - \frac{m-1}{\kappa} \end{aligned}
=⇒∂x∂[(1−rx)m−1exp(κx)]−r(m−1)(1−rx)m−2exp(κx)+(1−rx)m−1κexp(κx)=0−r(m−1)+(1−rx)κ=01−rx=κr(m−1)rx=1−κr(m−1)解: x0=r1−κm−1 因为
r
∈
(
0
,
1
)
,
x
∈
[
−
1
,
1
]
r \in (0,1), x \in [-1, 1]
r∈(0,1),x∈[−1,1], 所以
1
−
r
x
≥
0
1-rx \ge 0
1−rx≥0, 可知, 在经过简约之后,
∂
[
(
1
−
r
x
)
m
−
1
e
x
p
(
κ
x
)
]
∂
x
\frac{\partial [(1-rx)^{m-1} exp(\kappa x)]}{\partial x}
∂x∂[(1−rx)m−1exp(κx)] 的符号不变,
−
r
(
m
−
1
)
+
(
1
−
r
x
)
κ
-r(m-1) + (1-rx) \kappa
−r(m−1)+(1−rx)κ 是线性的, 单调递减. 那么
f
r
a
d
i
a
l
(
x
;
κ
,
m
)
e
(
x
,
b
)
\frac{f_{radial}(x;\kappa,m)}{e(x, b)}
e(x,b)fradial(x;κ,m) 先曾后减, 在
x
0
x_0
x0 处取得最大值.
代入之后, 求得 M M M, 为使接受率最大化(最小化 M M M), 对 r r r 求导: ∂ ( [ ( 1 − r x 0 ) m − 1 e x p ( κ x 0 ) ] ( 1 − r 2 ) m − 1 2 ) ∂ b = e x p ( − m − 1 κ ) ∂ ( [ r ( m − 1 ) κ ] m − 1 e x p ( κ r ) ( 1 − r 2 ) m − 1 2 ) ∂ b = ( m − 1 κ ) m − 1 e x p ( − m − 1 κ ) ∂ ( r m − 1 e x p ( κ r ) ( 1 − r 2 ) m − 1 2 ) ∂ b = 0 ⇒ = ∂ ( r m − 1 e x p ( κ r ) ( 1 − r 2 ) m − 1 2 ) ∂ b = 0 = m − 1 2 2 r ( 1 − r 2 ) 2 ( r 2 1 − r 2 ) m − 3 2 e x p ( κ r ) + ( r 2 1 − r 2 ) m − 1 2 − κ r 2 e x p ( κ r ) ⇒ m − 1 2 2 r ( 1 − r 2 ) 2 + r 2 1 − r 2 − κ r 2 = 0 r ( m − 1 ) ( 1 − r 2 ) 2 − κ 1 − r 2 = 0 r ( m − 1 ) − κ ( 1 − r 2 ) = 0 ⇒ κ r 2 + ( m − 1 ) r − κ = 0 解 : r 0 = − ( m − 1 ) + 4 κ 2 + ( m − 1 ) 2 2 κ ∈ ( 0 , 1 ) \begin{aligned} & \frac{\partial \left( \frac{[(1-rx_0)^{m-1} exp(\kappa x_0)]}{(1-r^2)^\frac{m-1}{2}} \right)}{\partial b} \\ =& exp(-\frac{m-1}{\kappa}) \frac{\partial \left( \frac{[\frac{r(m-1)}{\kappa}]^{m-1} exp(\frac{\kappa}{r})}{(1-r^2)^\frac{m-1}{2}} \right)}{\partial b} \\ =& \left(\frac{m-1}{\kappa}\right)^{m-1} exp(-\frac{m-1}{\kappa}) \frac{\partial \left( \frac{r^{m-1} exp(\frac{\kappa}{r})}{(1-r^2)^\frac{m-1}{2}} \right)}{\partial b} = 0 \\ \Rightarrow \\ =& \frac{\partial \left( \frac{r^{m-1} exp(\frac{\kappa}{r})}{(1-r^2)^\frac{m-1}{2}} \right)}{\partial b} = 0 \\ =& \frac{m-1}{2} \frac{2r}{(1-r^2)^2} (\frac{r^2}{1-r^2})^{\frac{m-3}{2}} exp(\frac{\kappa}{r}) + (\frac{r^2}{1-r^2})^{\frac{m-1}{2}} \frac{-\kappa}{r^2} exp(\frac{\kappa}{r}) \\ \Rightarrow \\ & \frac{m-1}{2} \frac{2r}{(1-r^2)^2} + \frac{r^2}{1-r^2} \frac{-\kappa}{r^2} = 0 \\ & \frac{r(m-1)}{(1-r^2)^2} - \frac{\kappa}{1-r^2} = 0 \\ & r(m-1) - \kappa(1-r^2) = 0 \\ \Rightarrow & \kappa r^2 + (m-1)r - \kappa = 0 \\ & 解:~ r_0 = \frac{-(m-1)+\sqrt{4\kappa^2+(m-1)^2}}{2\kappa} \in (0,1) \end{aligned} ==⇒==⇒⇒∂b∂((1−r2)2m−1[(1−rx0)m−1exp(κx0)])exp(−κm−1)∂b∂((1−r2)2m−1[κr(m−1)]m−1exp(rκ))(κm−1)m−1exp(−κm−1)∂b∂((1−r2)2m−1rm−1exp(rκ))=0∂b∂((1−r2)2m−1rm−1exp(rκ))=02m−1(1−r2)22r(1−r2r2)2m−3exp(rκ)+(1−r2r2)2m−1r2−κexp(rκ)2m−1(1−r2)22r+1−r2r2r2−κ=0(1−r2)2r(m−1)−1−r2κ=0r(m−1)−κ(1−r2)=0κr2+(m−1)r−κ=0解: r0=2κ−(m−1)+4κ2+(m−1)2∈(0,1)
好眼熟啊, 不就是上一个拒绝采样法的 x 0 x_0 x0 吗?
所有的简约都不改变 ∂ M ∂ r \frac{\partial M}{\partial r} ∂r∂M 的符号, 所以, M M M 随 r 先减后增, r 0 r_0 r0 是 M M M 的最小值点.
继续计算 x 0 x_0 x0: x 0 = 1 r 0 − m − 1 κ = 2 κ − ( m − 1 ) + 4 κ 2 + ( m − 1 ) 2 − m − 1 κ = 2 κ ( ( m − 1 ) + 4 κ 2 + ( m − 1 ) 2 ) ( − ( m − 1 ) + 4 κ 2 + ( m − 1 ) 2 ) ( ( m − 1 ) + 4 κ 2 + ( m − 1 ) 2 ) − m − 1 κ = 2 κ ( ( m − 1 ) + 4 κ 2 + ( m − 1 ) 2 ) − ( m − 1 ) 2 + 4 κ 2 + ( m − 1 ) 2 − m − 1 κ = ( m − 1 ) + 4 κ 2 + ( m − 1 ) 2 2 κ − 2 ( m − 1 ) 2 κ = − ( m − 1 ) + 4 κ 2 + ( m − 1 ) 2 2 κ = r 0 \begin{aligned} x_0 &= \frac{1}{r_0} - \frac{m-1}{\kappa} \\ &= \frac{2\kappa}{-(m-1)+\sqrt{4\kappa^2+(m-1)^2}} - \frac{m-1}{\kappa} \\ &= \frac{2\kappa ((m-1)+\sqrt{4\kappa^2+(m-1)^2})}{(-(m-1)+\sqrt{4\kappa^2+(m-1)^2})((m-1)+\sqrt{4\kappa^2+(m-1)^2})} - \frac{m-1}{\kappa} \\ &= \frac{2\kappa ((m-1)+\sqrt{4\kappa^2+(m-1)^2})}{-(m-1)^2+4\kappa^2+(m-1)^2} - \frac{m-1}{\kappa} \\ &= \frac{(m-1)+\sqrt{4\kappa^2+(m-1)^2}}{2\kappa} - \frac{2(m-1)}{2\kappa} \\ &= \frac{-(m-1)+\sqrt{4\kappa^2+(m-1)^2}}{2\kappa} \\ &= r_0 \end{aligned} x0=r01−κm−1=−(m−1)+4κ2+(m−1)22κ−κm−1=(−(m−1)+4κ2+(m−1)2)((m−1)+4κ2+(m−1)2)2κ((m−1)+4κ2+(m−1)2)−κm−1=−(m−1)2+4κ2+(m−1)22κ((m−1)+4κ2+(m−1)2)−κm−1=2κ(m−1)+4κ2+(m−1)2−2κ2(m−1)=2κ−(m−1)+4κ2+(m−1)2=r0 采样一个 u ∼ U n i f o r m ( 0 , 1 ) u \sim Uniform(0,1) u∼Uniform(0,1), y ∼ B e ( m − 1 2 , m − 1 2 ) ⇒ x = r 0 + ( 2 y − 1 ) 1 + r 0 ( 2 y − 1 ) ∼ e ( x , r 0 ) y \sim Be(\frac{m-1}{2},\frac{m-1}{2}) \Rightarrow x = \frac{r_0+(2y-1)}{1+r_0(2y-1)} \sim e(x,r_0) y∼Be(2m−1,2m−1)⇒x=1+r0(2y−1)r0+(2y−1)∼e(x,r0) 当 u < f r a d i a l ( x ; κ , m ) M ∗ e ( x , r 0 ) u \lt \frac{f_{radial}(x;\kappa,m)}{M*e(x,r_0)} u<M∗e(x,r0)fradial(x;κ,m) 时, 接受 x x x 为样本. f r a d i a l ( x ; κ , m ) M ∗ e ( x , b 0 ) = f r a d i a l ( x ; κ , m ) M f ( ( 1 − r 0 ) ( 1 + x ) 2 ( 1 − r 0 x ) ; m − 1 2 , m − 1 2 ) = ( κ / 2 ) ν Γ ( 1 2 ) Γ ( ν + 1 2 ) I ν ( κ ) 2 m − 2 B ( m − 1 2 , m − 1 2 ) ( 1 − r 0 2 ) m − 1 2 ( 1 − r 0 x ) m − 1 e x p ( κ x ) / M = ( 1 − r 0 x ) m − 1 e x p ( κ x ) ( 1 − r 0 x 0 ) m − 1 e x p ( κ x 0 ) = ( 1 − r 0 x 1 − r 0 x 0 ) m − 1 e x p ( κ ( x − x 0 ) ) l o g f r a d i a l ( x ; κ , m ) M ∗ e ( x , r 0 ) = ( m − 1 ) l o g ( 1 − r 0 x 1 − r 0 x 0 ) + κ ( x − x 0 ) = ( m − 1 ) l o g ( 1 − x 0 x ) − ( m − 1 ) l o g ( 1 − x 0 2 ) + κ ( x − x 0 ) = κ x + ( m − 1 ) l o g ( 1 − x 0 x ) − [ κ x 0 + ( m − 1 ) l o g ( 1 − x 0 2 ) ] = κ x + ( m − 1 ) l o g ( 1 − x 0 x ) − c \begin{aligned} \frac{f_{radial}(x;\kappa,m)}{M*e(x,b_0)} &= \frac{f_{radial}(x;\kappa,m)}{Mf(\frac{(1-r_0)(1+x)}{2(1-r_0x)}; \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{2^{m-2}\Beta(\frac{m-1}{2},\frac{m-1}{2})}{(1-r_0^2)^\frac{m-1}{2}} (1-r_0x)^{m-1} exp(\kappa x) / M \\ &= \frac{(1-r_0x)^{m-1} exp(\kappa x)}{(1-r_0x_0)^{m-1} exp(\kappa x_0)} \\ &= \left(\frac{1-r_0x}{1-r_0x_0}\right)^{m-1} exp(\kappa (x - x_0)) \\ log\frac{f_{radial}(x;\kappa,m)}{M*e(x,r_0)} &= (m-1) log\left(\frac{1-r_0x}{1-r_0x_0}\right) + \kappa (x - x_0) \\ &= (m-1)log(1-x_0x) - (m-1)log(1-x_0^2) + \kappa (x - x_0) \\ &= \kappa x + (m-1)log(1-x_0x) - [\kappa x_0 + (m-1)log(1-x_0^2)] \\ &= \kappa x + (m-1)log(1-x_0x) - c \end{aligned} M∗e(x,b0)fradial(x;κ,m)logM∗e(x,r0)fradial(x;κ,m)=Mf(2(1−r0x)(1−r0)(1+x);2m−1,2m−1)fradial(x;κ,m)=Γ(21)Γ(ν+21)Iν(κ)(κ/2)ν(1−r02)2m−12m−2B(2m−1,2m−1)(1−r0x)m−1exp(κx)/M=(1−r0x0)m−1exp(κx0)(1−r0x)m−1exp(κx)=(1−r0x01−r0x)m−1exp(κ(x−x0))=(m−1)log(1−r0x01−r0x)+κ(x−x0)=(m−1)log(1−x0x)−(m−1)log(1−x02)+κ(x−x0)=κx+(m−1)log(1−x0x)−[κx0+(m−1)log(1−x02)]=κx+(m−1)log(1−x0x)−c 即, 当 l o g u < κ x + ( m − 1 ) l o g ( 1 − x 0 x ) − c logu < \kappa x + (m-1)log(1-x_0x) - c logu<κx+(m−1)log(1−x0x)−c 时接受样本 x x x. 兜兜转转, 还是同一个判别式, 连 x 0 x_0 x0 都是一模一样的. 所以, 所有都一模一样的! 甚至: e ( x , b ) = 2 b m − 1 2 B ( m − 1 2 , m − 1 2 ) ( 1 − x 2 ) m − 3 2 [ ( 1 + b ) − ( 1 − b ) x ] m − 1 = 2 b m − 1 2 B ( m − 1 2 , m − 1 2 ) ( 1 − x 2 ) m − 3 2 ( 1 + b ) m − 1 ( 1 − 1 − b 1 + b x ) m − 1 令 r = 1 − b 1 + b ⇔ b = 1 − r 1 + r = 2 ( 4 b ( 1 + b ) 2 ) m − 1 2 2 ( m − 1 ) B ( m − 1 2 , m − 1 2 ) ( 1 − x 2 ) m − 3 2 ( 1 − r x ) m − 1 = ( ( 1 + b ) 2 − ( 1 − b ) 2 ( 1 + b ) 2 ) m − 1 2 2 m − 2 B ( m − 1 2 , m − 1 2 ) ( 1 − x 2 ) m − 3 2 ( 1 − r x ) m − 1 = ( 1 − r 2 ) m − 1 2 2 m − 2 B ( m − 1 2 , m − 1 2 ) ( 1 − x 2 ) m − 3 2 ( 1 − r x ) m − 1 \begin{aligned} e(x,b) =& \frac{2b^{\frac{m-1}{2}}}{\Beta(\frac{m-1}{2},\frac{m-1}{2})} \frac{(1-x^2)^{\frac{m-3}{2}}}{[(1+b)-(1-b)x]^{m-1}} \\ =& \frac{2b^{\frac{m-1}{2}}}{\Beta(\frac{m-1}{2},\frac{m-1}{2})} \frac{(1-x^2)^{\frac{m-3}{2}}}{(1+b)^{m-1}\left(1-\frac{1-b}{1+b}x\right)^{m-1}} \\ 令 ~ r=& \frac{1-b}{1+b} \Leftrightarrow b= \frac{1-r}{1+r} \\ =& \frac{ 2\left(\frac{4b}{(1+b)^2}\right)^\frac{m-1}{2} }{2^{(m-1)}\Beta(\frac{m-1}{2},\frac{m-1}{2})} \frac{(1-x^2)^{\frac{m-3}{2}}}{\left(1-rx\right)^{m-1}} \\ =& \frac{ \left(\frac{(1+b)^2 - (1-b)^2}{(1+b)^2}\right)^\frac{m-1}{2} }{2^{m-2}\Beta(\frac{m-1}{2},\frac{m-1}{2})} \frac{(1-x^2)^{\frac{m-3}{2}}}{\left(1-rx\right)^{m-1}} \\ =& \frac{(1-r^2)^\frac{m-1}{2}}{2^{m-2}\Beta(\frac{m-1}{2},\frac{m-1}{2})} \frac{(1-x^2)^{\frac{m-3}{2}}}{\left(1-rx\right)^{m-1}} \end{aligned} e(x,b)==令 r====B(2m−1,2m−1)2b2m−1[(1+b)−(1−b)x]m−1(1−x2)2m−3B(2m−1,2m−1)2b2m−1(1+b)m−1(1−1+b1−bx)m−1(1−x2)2m−31+b1−b⇔b=1+r1−r2(m−1)B(2m−1,2m−1)2((1+b)24b)2m−1(1−rx)m−1(1−x2)2m−32m−2B(2m−1,2m−1)((1+b)2(1+b)2−(1−b)2)2m−1(1−rx)m−1(1−x2)2m−32m−2B(2m−1,2m−1)(1−r2)2m−1(1−rx)m−1(1−x2)2m−3 建议概率密度函数是一模一样的.
但是, 当将
r
=
1
−
b
1
+
b
r = \frac{1-b}{1+b}
r=1+b1−b 代入
x
x
x 时, 计算:
x
=
r
0
+
(
2
y
−
1
)
1
+
r
0
(
2
y
−
1
)
=
x
0
+
(
2
y
−
1
)
1
+
x
0
(
2
y
−
1
)
=
1
−
b
0
1
+
b
0
+
(
2
y
−
1
)
1
+
1
−
b
0
1
+
b
0
(
2
y
−
1
)
=
1
−
b
0
+
(
1
+
b
0
)
(
2
y
−
1
)
1
+
b
0
+
(
1
−
b
0
)
(
2
y
−
1
)
=
1
−
b
0
+
2
y
+
2
b
0
y
−
1
−
b
0
1
+
b
0
+
2
y
−
2
b
0
y
−
1
+
b
0
=
b
0
−
(
1
+
b
0
)
y
−
b
0
−
(
1
−
b
0
)
y
\begin{aligned} x &= \frac{r_0+(2y-1)}{1+r_0(2y-1)} \\ &= \frac{x_0+(2y-1)}{1+x_0(2y-1)} \\ &= \frac{\frac{1-b_0}{1+b_0}+(2y-1)}{1+\frac{1-b_0}{1+b_0}(2y-1)} \\ &= \frac{1-b_0+(1+b_0)(2y-1)}{1+b_0+(1-b_0)(2y-1)} \\ &= \frac{1 - b_0 + 2y + 2b_0y - 1 - b_0}{1 + b_0 + 2y - 2b_0y - 1 + b_0} \\ &= \frac{b_0 - (1+b_0)y}{-b_0 - (1 - b_0)y} \\ \end{aligned}
x=1+r0(2y−1)r0+(2y−1)=1+x0(2y−1)x0+(2y−1)=1+1+b01−b0(2y−1)1+b01−b0+(2y−1)=1+b0+(1−b0)(2y−1)1−b0+(1+b0)(2y−1)=1+b0+2y−2b0y−1+b01−b0+2y+2b0y−1−b0=−b0−(1−b0)yb0−(1+b0)y 怎么也化不成
1
−
(
1
+
b
0
)
y
1
−
(
1
−
b
0
)
y
\frac{1 - (1+b_0)y}{1 - (1 - b_0)y}
1−(1−b0)y1−(1+b0)y, 画出图像:
才发现, 两者不是相等的, 而是关于
y
=
0.5
y=0.5
y=0.5 对称. 不过话说回来,
B
e
(
y
,
m
−
1
2
,
m
−
1
2
)
Be(y, \frac{m-1}{2}, \frac{m-1}{2})
Be(y,2m−1,2m−1) 分布确实是关于
y
=
0.5
y=0.5
y=0.5 对称的. 如果把
y
y
y 换成
1
−
y
1-y
1−y 继续:
b
0
−
(
1
+
b
0
)
(
1
−
y
)
−
b
0
−
(
1
−
b
0
)
(
1
−
y
)
=
b
0
−
(
1
+
b
0
)
+
(
1
+
b
0
)
y
−
b
0
−
(
1
−
b
0
)
+
(
1
−
b
0
)
y
=
−
1
+
(
1
+
b
0
)
y
−
1
+
(
1
−
b
0
)
y
=
1
−
(
1
+
b
0
)
y
1
−
(
1
−
b
0
)
y
\begin{aligned} \frac{b_0 - (1+b_0)(1-y)}{-b_0 - (1 - b_0)(1-y)} &= \frac{b_0 - (1+b_0) + (1+b_0)y}{-b_0 - (1 - b_0) + (1 - b_0)y} \\ &= \frac{-1 + (1+b_0)y}{-1 + (1 - b_0)y} \\ &= \frac{1 - (1+b_0)y}{1 - (1 - b_0)y} \end{aligned}
−b0−(1−b0)(1−y)b0−(1+b0)(1−y)=−b0−(1−b0)+(1−b0)yb0−(1+b0)+(1+b0)y=−1+(1−b0)y−1+(1+b0)y=1−(1−b0)y1−(1+b0)y OK 了.
参考文献
[3] Fast Python Sampler of the von Mises Fisher Distribution.
![[管理者与领导者-174] ：人际网络-1- 网络概述，是由一个个人组成的网络，每个节点是“人”](https://img-blog.csdnimg.cn/direct/12aacd25422249ceb949228373bcc02e.png)
