2. Relation to Normal Distribution
疑问:有没有不各向同性的 vMF?
答:应该是没有的,如果想让各方向偏离中心的速度不一致,则协方差矩阵不为
I
\bm{I}
I 的倍数. 正态分布的概率密度函数为:
f
(
x
)
=
1
(
2
π
)
p
∣
Σ
∣
e
−
1
2
(
x
−
μ
)
⊺
Σ
−
1
(
x
−
μ
)
\begin{aligned} f(\bm{x}) = \frac{1}{\sqrt{(2\pi)^p |\Sigma|}} e^{-\frac{1}{2} (\bm{x}-\bm{\mu})^\intercal\Sigma^{-1}(\bm{x}-\bm{\mu})} \end{aligned}
f(x)=(2π)p∣Σ∣1e−21(x−μ)⊺Σ−1(x−μ) 类比上面的推导,我们需要得出形似:
G
p
(
x
;
μ
,
κ
)
=
C
(
p
,
κ
,
r
)
e
x
p
(
κ
r
μ
⊺
r
Σ
−
1
x
)
\begin{aligned} G_p(\bm{x}; \bm{\mu}, \kappa) &= C(p,\kappa,r) exp\left(\kappa r \frac{\bm{\mu}^\intercal}{r} \Sigma^{-1} \bm{x} \right) \end{aligned}
Gp(x;μ,κ)=C(p,κ,r)exp(κrrμ⊺Σ−1x) 的东西,所以,必要的,需要:
x
⊺
Σ
−
1
x
=
c
o
n
s
t
\begin{aligned} \bm{x}^{\intercal} \Sigma^{-1} \bm{x} &= const \end{aligned}
x⊺Σ−1x=const 我们知道,要想让
x
⊺
Σ
−
1
x
=
c
o
n
s
t
\bm{x}^{\intercal} \Sigma^{-1} \bm{x} = const
x⊺Σ−1x=const 表示为球,则必须使
Σ
−
1
=
α
I
\Sigma^{-1}=\alpha \bm{I}
Σ−1=αI,所以,假设不了“协方差矩阵不为
I
\bm{I}
I 的倍”, 也就不可能存在不各向同性的 vMF.
[注]: x ⊺ Σ − 1 x = c o n s t \bm{x}^{\intercal} \Sigma^{-1} \bm{x} = const x⊺Σ−1x=const,表示一般的超椭球。
等等!对于
x
\bm{x}
x 来说, 即使它不是单位向量,也代表了一个方向,
x
∥
x
∥
\frac{\bm{x}}{\|\bm{x}\|}
∥x∥x 的分布会是非各向同性的 vMF 吗?有可能是哎!
假设
y
=
[
y
1
y
2
]
=
[
1
0
0
2
]
[
x
1
x
2
]
=
A
x
\bm{y} = \begin{bmatrix} y_1 \\ y_2 \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 0 & 2 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \bm{Ax}
y=[y1y2]=[1002][x1x2]=Ax,
y
\bm{y}
y 代表单位圆的话,则
x
\bm{x}
x 在椭圆
x
1
2
+
(
2
x
2
)
2
=
1
x_1^2 + (2x_2)^2 = 1
x12+(2x2)2=1 上,如下图:
可以看到,B 和 D 点是对应的,那么弧 AB 和 AD 的“点数”应该是一样的,而弧 AB 对应的方向却在弧 AC 上,即 AD 段对应的方向,被压缩在了 AC 段,如果采样
y
\bm{y}
y 的话,对应的
x
\bm{x}
x 方向会更集中趋于椭圆的长轴,实现了“非各向同性”,只不过分布在椭圆上,且采样后需要归一化处理。
5. 采样
t t t 的概率密度函数推导
直接将 t = μ ⊺ x t = \bm{\mu}^\intercal\bm{x} t=μ⊺x 代入 f p ( x ; μ , κ ) f_p(\bm{x}; \bm{\mu}, \kappa) fp(x;μ,κ), 得: f p ( x ; μ , κ ) = C p ( κ ) e x p ( κ μ ⊺ x ) = C p ( κ ) e x p ( κ t ) t ∈ [ − 1 , 1 ] = C p ( κ ) e x p ( κ c o s θ ) θ ∈ [ 0 , π ] \begin{aligned} f_p(\bm{x}; \bm{\mu}, \kappa) &= C_p(\kappa) exp(\kappa \bm{\mu}^\intercal\bm{x}) & \\ &= C_p(\kappa) exp(\kappa t) & t \in [-1, 1] \\ &= C_p(\kappa) exp(\kappa cos\theta) & \theta \in [0, \pi] \end{aligned} fp(x;μ,κ)=Cp(κ)exp(κμ⊺x)=Cp(κ)exp(κt)=Cp(κ)exp(κcosθ)t∈[−1,1]θ∈[0,π] 注意这是 x \bm{x} x 一点处的概率密度. 沿着 t t t 处的切子球求积分, 以得到 t t t 或 θ \theta θ 处的整个概率密度: ∫ 切子球 f p ( x ; μ , κ ) d s = ∫ 切子球 C p ( κ ) e x p ( κ x ) d s = C p ( κ ) e x p ( κ t ) 2 π p − 1 2 Γ ( p − 1 2 ) ( 1 − t 2 ) p − 2 2 S p − 2 的表面积 ∝ r p − 2 = C p ( κ ) e x p ( κ c o s θ ) 2 π p − 1 2 Γ ( p − 1 2 ) s i n p − 2 θ \begin{aligned} \int_{切子球} f_p(\bm{x}; \bm{\mu}, \kappa) ds &= \int_{切子球} C_p(\kappa) exp(\kappa \bm{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})} sin^{p-2}\theta & \end{aligned} ∫切子球fp(x;μ,κ)ds=∫切子球Cp(κ)exp(κx)ds=Cp(κ)exp(κt)Γ(2p−1)2π2p−1(1−t2)2p−2=Cp(κ)exp(κcosθ)Γ(2p−1)2π2p−1sinp−2θSp−2的表面积∝rp−2 根据 n-sphere - Wikipedia, 切子球 S p − 2 S^{p-2} Sp−2 的表面积 S p − 2 = 2 π p − 1 2 Γ ( p − 1 2 ) r p − 2 S_{p-2} = \frac{2\pi^{\frac{p-1}{2}}}{\Gamma(\frac{p-1}{2})} r^{p-2} Sp−2=Γ(2p−1)2π2p−1rp−2, 再沿 t t t 或 θ \theta θ 积分: ∫ 0 π C p ( κ ) e x p ( κ c o s θ ) 2 π p − 1 2 Γ ( p − 1 2 ) s i n p − 2 θ d θ = C p ( κ ) 2 π p − 1 2 Γ ( p − 1 2 ) ∫ 1 − 1 e x p ( κ t ) ( 1 − t 2 ) p − 2 2 ( − 1 1 − t 2 d t ) ∵ c o s 0 = 1 , c o s π = − 1 = C p ( κ ) 2 π p − 1 2 Γ ( p − 1 2 ) ∫ − 1 1 e x p ( κ t ) ( 1 − t 2 ) p − 3 2 d t = C p ( κ ) 2 π ν + 1 2 Γ ( ν + 1 2 ) ∫ − 1 1 e x p ( κ t ) ( 1 − t 2 ) ν − 1 2 d t 令 ν = p 2 − 1 \begin{aligned} & \int_0^\pi C_p(\kappa) exp(\kappa cos\theta) \frac{2\pi^{\frac{p-1}{2}}}{\Gamma(\frac{p-1}{2})} sin^{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{aligned} ===∫0πCp(κ)exp(κcosθ)Γ(2p−1)2π2p−1sinp−2θdθCp(κ)Γ(2p−1)2π2p−1∫1−1exp(κt)(1−t2)2p−2(1−t2−1dt)Cp(κ)Γ(2p−1)2π2p−1∫−11exp(κt)(1−t2)2p−3dtCp(κ)Γ(ν+21)2πν+21∫−11exp(κt)(1−t2)ν−21dt∵cos0=1, cosπ=−1令 ν=2p−1 那么, 将 ∫ − 1 1 e x p ( κ t ) ( 1 − t 2 ) ν − 1 2 d t = Γ ( 1 2 ) Γ ( ν + 1 2 ) I ν ( κ ) ( κ / 2 ) ν [ I ν 的公式 ] = Γ ( 1 2 ) Γ ( p − 1 2 ) I p 2 − 1 ( κ ) ( κ / 2 ) p 2 − 1 \begin{aligned} \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{aligned} ∫−11exp(κt)(1−t2)ν−21dt=(κ/2)νΓ(21)Γ(ν+21)Iν(κ)=(κ/2)2p−1Γ(21)Γ(2p−1)I2p−1(κ)[Iν的公式] 代入, 得: C p ( κ ) 2 π ν + 1 2 Γ ( ν + 1 2 ) ∫ − 1 1 e x p ( κ t ) ( 1 − t 2 ) ν − 1 2 d t = C p ( κ ) 2 π ν + 1 2 Γ ( ν + 1 2 ) Γ ( 1 2 ) Γ ( ν + 1 2 ) I ν ( κ ) ( κ / 2 ) ν = 2 π ν + 1 I ν ( κ ) ( κ / 2 ) ν C p ( κ ) Γ ( 1 2 ) = π 1 2 = 2 π p 2 I p 2 − 1 ( κ ) ( κ / 2 ) p 2 − 1 C p ( κ ) = ( 2 π ) p 2 I p 2 − 1 ( κ ) κ p 2 − 1 C p ( κ ) = 1 ( 2 ) ⟺ C p ( κ ) = κ p 2 − 1 ( 2 π ) p 2 I p 2 − 1 ( κ ) [ 与 W i k i p e d i a 一致 ] \begin{aligned} & 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) \\ \Longleftrightarrow & \\ C_p(\kappa) =& \frac{\kappa^{\frac{p}{2}-1}}{(2\pi)^{\frac{p}{2}} I_{\frac{p}{2}-1}(\kappa)} & [与 Wikipedia 一致] \end{aligned} ====⟺Cp(κ)=Cp(κ)Γ(ν+21)2πν+21∫−11exp(κt)(1−t2)ν−21dtCp(κ)Γ(ν+21)2πν+21(κ/2)νΓ(21)Γ(ν+21)Iν(κ)(κ/2)ν2πν+1Iν(κ)Cp(κ)(κ/2)2p−12π2pI2p−1(κ)Cp(κ)κ2p−1(2π)2pI2p−1(κ)Cp(κ)=1(2π)2pI2p−1(κ)κ2p−1Γ(21)=π21(2)[与Wikipedia一致] 所以, 由式 ( 2 ) (2) (2) 可得 t t t 的概率密度函数为: f r a d i a l ( t ; κ , p ) = C p ( κ ) 2 π ν + 1 2 Γ ( ν + 1 2 ) e x p ( κ t ) ( 1 − t 2 ) ν − 1 2 ν = p 2 − 1 = ( κ / 2 ) ν Γ ( 1 2 ) Γ ( ν + 1 2 ) I ν ( κ ) e x p ( κ t ) ( 1 − t 2 ) ν − 1 2 [ 与 W i k i p e d i a 一致 ] = ( κ / 2 ) p 2 − 1 π Γ ( p − 1 2 ) I p 2 − 1 ( κ ) e x p ( κ t ) ( 1 − t 2 ) p − 3 2 \begin{aligned} 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{aligned} fradial(t;κ,p)=Cp(κ)Γ(ν+21)2πν+21exp(κt)(1−t2)ν−21=Γ(21)Γ(ν+21)Iν(κ)(κ/2)νexp(κt)(1−t2)ν−21=πΓ(2p−1)I2p−1(κ)(κ/2)2p−1exp(κt)(1−t2)2p−3ν=2p−1[与Wikipedia一致]
单位球上的均匀分布采样
实际上, 从标准正态分布 N o r m a l ( 0 , 1 ) Normal(0,1) Normal(0,1) 中采样 p p p 个数, 组成向量 v \bm{v} v, 归一化, 就实现了球 S p − 1 S_{p-1} Sp−1 上的均匀采样, 并不需要多元正态分布 M u l t i v a r i a t e N o r m a l ( 0 , I ( p ) ) MultivariateNormal(0, \bm{I}^{(p)}) MultivariateNormal(0,I(p)). 因为 p p p 个独立采样自 N o r m a l ( 0 , 1 ) Normal(0,1) Normal(0,1) 的 p p p 个数概率是: ∏ i = 1 p 1 2 π e x p ( − x i 2 2 ) = 1 ( 2 π ) p 2 e x p ( − ∑ i = 1 p x i 2 2 ) \begin{aligned} & \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{aligned} =i=1∏p2π1exp(−2xi2)(2π)2p1exp(−2∑i=1pxi2) 恰好是 p p p 元标准正态分布的概率密度.
5.2 拒绝采样
为了与原论文一致, 本小节使用 m m m 代替 p p p.
首先, 设
Y
∼
B
e
(
α
,
β
)
Y \sim Be(\alpha, \beta)
Y∼Be(α,β), 概率密度函数为:
f
(
y
;
α
,
β
)
=
y
α
−
1
(
1
−
y
)
β
−
1
∫
0
1
u
α
−
1
(
1
−
u
)
β
−
1
d
u
=
Γ
(
α
+
β
)
Γ
(
α
)
Γ
(
β
)
y
α
−
1
(
1
−
y
)
β
−
1
=
1
B
(
α
,
β
)
y
α
−
1
(
1
−
y
)
β
−
1
\begin{aligned} 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}{\Beta(\alpha,\beta)} y^{\alpha-1}(1-y)^{\beta-1} \end{aligned}
f(y;α,β)=∫01uα−1(1−u)β−1duyα−1(1−y)β−1=Γ(α)Γ(β)Γ(α+β)yα−1(1−y)β−1=B(α,β)1yα−1(1−y)β−1 其中
y
∈
(
0
,
1
)
,
α
,
β
>
0
y \in (0,1), \alpha, \beta > 0
y∈(0,1),α,β>0. 令
X
=
1
−
(
1
+
b
)
Y
1
−
(
1
−
b
)
Y
⇔
Y
=
1
−
X
(
1
+
b
)
−
(
1
−
b
)
X
b
∈
(
0
,
1
)
\begin{aligned} X = \frac{1-(1+b)Y}{1-(1-b)Y} &\Leftrightarrow Y = \frac{1-X}{(1+b)-(1-b)X} \\ b &\in (0,1) \end{aligned}
X=1−(1−b)Y1−(1+b)Yb⇔Y=(1+b)−(1−b)X1−X∈(0,1) 代入上述概率密度函数
f
(
y
;
α
,
β
)
f(y; \alpha, \beta)
f(y;α,β), 并令
α
=
β
=
m
−
1
2
\alpha=\beta=\frac{m-1}{2}
α=β=2m−1 得:
f
(
1
−
x
(
1
+
b
)
−
(
1
−
b
)
x
;
m
−
1
2
,
m
−
1
2
)
=
1
B
(
m
−
1
2
,
m
−
1
2
)
(
(
1
−
x
(
1
+
b
)
−
(
1
−
b
)
x
)
∗
(
1
−
1
−
x
(
1
+
b
)
−
(
1
−
b
)
x
)
)
m
−
1
2
−
1
=
1
B
(
m
−
1
2
,
m
−
1
2
)
(
(
1
−
x
(
1
+
b
)
−
(
1
−
b
)
x
)
∗
(
(
1
+
b
)
−
(
1
−
b
)
x
−
(
1
−
x
)
(
1
+
b
)
−
(
1
−
b
)
x
)
)
m
−
3
2
=
1
B
(
m
−
1
2
,
m
−
1
2
)
(
(
1
−
x
(
1
+
b
)
−
(
1
−
b
)
x
)
∗
(
b
(
1
+
x
)
(
1
+
b
)
−
(
1
−
b
)
x
)
)
m
−
3
2
=
1
B
(
m
−
1
2
,
m
−
1
2
)
(
b
(
1
−
x
2
)
[
(
1
+
b
)
−
(
1
−
b
)
x
]
2
)
m
−
3
2
=
b
m
−
3
2
B
(
m
−
1
2
,
m
−
1
2
)
(
1
−
x
2
)
m
−
3
2
[
(
1
+
b
)
−
(
1
−
b
)
x
]
m
−
3
=
e
(
x
,
b
)
其中
x
∈
(
−
1
,
1
)
,
b
∈
(
0
,
1
)
\begin{aligned} &f(\frac{1-x}{(1+b)-(1-b)x}; \frac{m-1}{2}, \frac{m-1}{2}) \\ =& \frac{1}{\Beta(\frac{m-1}{2},\frac{m-1}{2})} \left(\left(\frac{1-x}{(1+b)-(1-b)x}\right) * \left(1-\frac{1-x}{(1+b)-(1-b)x}\right)\right)^{\frac{m-1}{2}-1} \\ =& \frac{1}{\Beta(\frac{m-1}{2},\frac{m-1}{2})} \left(\left(\frac{1-x}{(1+b)-(1-b)x}\right) * \left(\frac{(1+b)-(1-b)x - (1-x)}{(1+b)-(1-b)x}\right)\right)^{\frac{m-3}{2}} \\ =& \frac{1}{\Beta(\frac{m-1}{2},\frac{m-1}{2})} \left(\left(\frac{1-x}{(1+b)-(1-b)x}\right) * \left(\frac{b(1+x)}{(1+b)-(1-b)x}\right)\right)^{\frac{m-3}{2}} \\ =& \frac{1}{\Beta(\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}}}{\Beta(\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{aligned}
======f((1+b)−(1−b)x1−x;2m−1,2m−1)B(2m−1,2m−1)1(((1+b)−(1−b)x1−x)∗(1−(1+b)−(1−b)x1−x))2m−1−1B(2m−1,2m−1)1(((1+b)−(1−b)x1−x)∗((1+b)−(1−b)x(1+b)−(1−b)x−(1−x)))2m−3B(2m−1,2m−1)1(((1+b)−(1−b)x1−x)∗((1+b)−(1−b)xb(1+x)))2m−3B(2m−1,2m−1)1([(1+b)−(1−b)x]2b(1−x2))2m−3B(2m−1,2m−1)b2m−3[(1+b)−(1−b)x]m−3(1−x2)2m−3e(x,b) 其中 x∈(−1,1),b∈(0,1) 与论文 Computer Generation of Distributions on the m-spher [1] 不一致, 但我不知道问题出在哪里. 论文里明明说:
却给出了:
由拒绝采样法的
f
r
a
d
i
a
l
(
x
;
κ
,
m
)
≤
M
∗
e
(
x
,
b
)
f_{radial}(x;\kappa,m) \le M*e(x,b)
fradial(x;κ,m)≤M∗e(x,b), 计算:
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
)
b
m
−
3
2
B
(
m
−
1
2
,
m
−
1
2
)
(
1
−
x
2
)
m
−
3
2
[
(
1
+
b
)
−
(
1
−
b
)
x
]
m
−
3
=
max
x
(
κ
/
2
)
ν
Γ
(
1
2
)
Γ
(
ν
+
1
2
)
I
ν
(
κ
)
e
x
p
(
κ
x
)
b
m
−
3
2
B
(
m
−
1
2
,
m
−
1
2
)
1
[
(
1
+
b
)
−
(
1
−
b
)
x
]
m
−
3
=
max
x
(
κ
/
2
)
ν
Γ
(
1
2
)
Γ
(
ν
+
1
2
)
I
ν
(
κ
)
B
(
m
−
1
2
,
m
−
1
2
)
b
m
−
3
2
[
(
1
+
b
)
−
(
1
−
b
)
x
]
m
−
3
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{b^{\frac{m-3}{2}}}{\Beta(\frac{m-1}{2},\frac{m-1}{2})} \frac{(1-x^2)^{\frac{m-3}{2}}}{[(1+b)-(1-b)x]^{m-3}} } \\ &= \max_x \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}}}{\Beta(\frac{m-1}{2},\frac{m-1}{2})} \frac{1}{[(1+b)-(1-b)x]^{m-3}} } \\ &= \max_x \frac{(\kappa/2)^\nu}{\Gamma({\frac{1}{2}})\Gamma(\nu+{\frac{1}{2}})I_{\nu}(\kappa)} \frac{\Beta(\frac{m-1}{2},\frac{m-1}{2})}{b^{\frac{m-3}{2}}} [(1+b)-(1-b)x]^{m-3} exp(\kappa x) \end{aligned}
M=xmaxe(x,b)fradial(x;κ,m)=xmaxB(2m−1,2m−1)b2m−3[(1+b)−(1−b)x]m−3(1−x2)2m−3Γ(21)Γ(ν+21)Iν(κ)(κ/2)ν(1−x2)2m−3exp(κx)=xmaxB(2m−1,2m−1)b2m−3[(1+b)−(1−b)x]m−31Γ(21)Γ(ν+21)Iν(κ)(κ/2)νexp(κx)=xmaxΓ(21)Γ(ν+21)Iν(κ)(κ/2)νb2m−3B(2m−1,2m−1)[(1+b)−(1−b)x]m−3exp(κx) 接下来求极值点:
∂
(
[
(
1
+
b
)
−
(
1
−
b
)
x
]
m
−
3
e
x
p
(
κ
x
)
)
∂
x
=
(
m
−
3
)
[
(
1
+
b
)
−
(
1
−
b
)
x
]
m
−
4
∗
[
−
(
1
−
b
)
]
e
x
p
(
κ
x
)
+
[
(
1
+
b
)
−
(
1
−
b
)
x
]
m
−
3
κ
e
x
p
(
κ
x
)
=
0
⇒
(
m
−
3
)
∗
[
−
(
1
−
b
)
]
+
[
(
1
+
b
)
−
(
1
−
b
)
x
]
κ
=
0
(
1
+
b
)
−
(
1
−
b
)
x
=
(
m
−
3
)
(
1
−
b
)
κ
(
1
−
b
)
x
=
(
1
+
b
)
−
(
m
−
3
)
(
1
−
b
)
κ
解
:
x
0
=
1
+
b
1
−
b
−
m
−
3
κ
\begin{aligned} & \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} * [-(1-b)] exp(\kappa x) \\ +& [(1+b)-(1-b)x]^{m-3} \kappa exp(\kappa x) = 0 \\ \Rightarrow \\ & (m-3) * [-(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{aligned}
=+⇒∂x∂([(1+b)−(1−b)x]m−3exp(κx))(m−3)[(1+b)−(1−b)x]m−4∗[−(1−b)]exp(κx)[(1+b)−(1−b)x]m−3κexp(κx)=0(m−3)∗[−(1−b)]+[(1+b)−(1−b)x]κ=0(1+b)−(1−b)x=κ(m−3)(1−b)(1−b)x=(1+b)−κ(m−3)(1−b)解: x0=1−b1+b−κm−3 因为
b
∈
(
0
,
1
)
,
x
∈
[
−
1
,
1
]
b \in (0,1), x \in [-1, 1]
b∈(0,1),x∈[−1,1], 所以
(
1
+
b
)
−
(
1
−
b
)
x
≥
0
(1+b)-(1-b)x \ge 0
(1+b)−(1−b)x≥0, 可知, 在经过简约之后,
∂
(
[
(
1
+
b
)
−
(
1
−
b
)
x
]
m
−
3
e
x
p
(
κ
x
)
)
∂
x
\frac{\partial \left([(1+b)-(1-b)x]^{m-3} exp(\kappa x)\right)}{\partial x}
∂x∂([(1+b)−(1−b)x]m−3exp(κx)) 的符号不变,
(
m
−
3
)
∗
[
−
(
1
−
b
)
]
+
[
(
1
+
b
)
−
(
1
−
b
)
x
]
κ
(m-3) * [-(1-b)] + [(1+b)-(1-b)x] \kappa
(m−3)∗[−(1−b)]+[(1+b)−(1−b)x]κ 是线性的, 单调递减. 那么
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), 对 b b b 求导: ∂ ( [ ( 1 + b ) − ( 1 − b ) x 0 ] m − 3 e x p ( κ x 0 ) b m − 3 2 ) ∂ b = e x p ( − m − 3 κ ) ∂ ( [ ( m − 3 ) ( 1 − b ) κ ] m − 3 e x p ( κ 1 + b 1 − b ) b m − 3 2 ) ∂ b = ( m − 3 κ ) m − 3 e x p ( − m − 3 κ ) ∂ ( ( 1 − b ) m − 3 e x p ( κ 1 + b 1 − b ) b m − 3 2 ) ∂ b = 0 ⇒ = ∂ ( ( 1 − b ) m − 3 e x p ( κ 1 + b 1 − b ) b m − 3 2 ) ∂ b = 0 = m − 3 2 ( 1 − 1 b 2 ) ( b + 1 b − 2 ) m − 5 2 e x p ( κ 1 + b 1 − b ) + ( b + 1 b − 2 ) m − 3 2 2 κ ( 1 − b ) 2 e x p ( κ 1 + b 1 − b ) ⇒ m − 3 2 ( 1 − 1 b 2 ) + ( b + 1 b − 2 ) 2 κ ( 1 − b ) 2 = 0 ( m − 3 ) ( 1 − 1 b 2 ) + 4 κ b = 0 ( m − 3 ) ( b 2 − 1 ) + 4 κ b = 0 = ( m − 3 ) b 2 + 4 κ b − ( m − 3 ) 解 : b 0 = − 2 κ + 4 κ 2 + ( m − 3 ) 2 m − 3 ∈ ( 0 , 1 ) \begin{aligned} & \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{aligned} ==⇒==⇒=∂b∂(b2m−3[(1+b)−(1−b)x0]m−3exp(κx0))exp(−κm−3)∂b∂(b2m−3[κ(m−3)(1−b)]m−3exp(κ1−b1+b))(κm−3)m−3exp(−κm−3)∂b∂(b2m−3(1−b)m−3exp(κ1−b1+b))=0∂b∂(b2m−3(1−b)m−3exp(κ1−b1+b))=02m−3(1−b21)(b+b1−2)2m−5exp(κ1−b1+b)+(b+b1−2)2m−3(1−b)22κexp(κ1−b1+b)2m−3(1−b21)+(b+b1−2)(1−b)22κ=0(m−3)(1−b21)+b4κ=0(m−3)(b2−1)+4κb=0(m−3)b2+4κb−(m−3)解: b0=m−3−2κ+4κ2+(m−3)2∈(0,1)
如果把 ( m − 3 ) > 0 (m-3) \gt 0 (m−3)>0 和 2 κ > 0 2\kappa \gt 0 2κ>0 都看作三角形的直角边, 4 κ 2 + ( m − 3 ) 2 \sqrt{4\kappa^2+(m-3)^2} 4κ2+(m−3)2 为斜边, 则 b 0 = s e c θ − t a n θ b_0 = sec\theta - tan\theta b0=secθ−tanθ, 又 θ ∈ ( 0 , π 2 ) \theta \in (0, \frac{\pi}{2}) θ∈(0,2π), 故 b 0 ∈ ( 0 , 1 ) b_0 \in (0,1) b0∈(0,1)
所有的简约都不改变 ∂ M ∂ b \frac{\partial M}{\partial b} ∂b∂M 的符号, 所以, M M M 随 b 先减后增, b 0 b_0 b0 是 M M M 的最小值点.
以上所有的计算中, 对于原论文中的
e
(
x
,
b
)
e(x,b)
e(x,b) 和本博文计算的
e
(
x
,
b
)
e(x,b)
e(x,b), 差别只在于
(
m
−
1
)
(m-1)
(m−1) 和
(
m
−
3
)
(m-3)
(m−3), 其他都一样. 继续计算
x
0
x_0
x0:
1
+
b
0
=
1
+
−
2
κ
+
4
κ
2
+
(
m
−
3
)
2
m
−
3
=
m
−
3
−
2
κ
+
4
κ
2
+
(
m
−
3
)
2
m
−
3
=
m
−
3
+
4
κ
2
+
(
m
−
3
)
2
−
2
κ
m
−
3
1
−
b
0
=
1
−
−
2
κ
+
4
κ
2
+
(
m
−
3
)
2
m
−
3
=
m
−
3
+
2
κ
−
4
κ
2
+
(
m
−
3
)
2
m
−
3
=
m
−
3
+
2
κ
−
4
κ
2
+
(
m
−
3
)
2
m
−
3
x
0
=
1
+
b
0
1
−
b
0
−
m
−
3
κ
=
m
−
3
−
2
κ
+
4
κ
2
+
(
m
−
3
)
2
m
−
3
+
2
κ
−
4
κ
2
+
(
m
−
3
)
2
−
m
−
3
κ
=
(
m
−
3
+
4
κ
2
+
(
m
−
3
)
2
−
2
κ
)
(
m
−
3
+
4
κ
2
+
(
m
−
3
)
2
+
2
κ
)
(
m
−
3
+
2
κ
)
2
−
(
4
κ
2
+
(
m
−
3
)
2
)
−
m
−
3
κ
=
(
m
−
3
+
4
κ
2
+
(
m
−
3
)
2
)
2
−
4
κ
2
(
m
−
3
)
2
+
4
κ
2
+
4
κ
(
m
−
3
)
−
(
4
κ
2
+
(
m
−
3
)
2
)
−
m
−
3
κ
=
(
m
−
3
)
2
+
4
κ
2
+
(
m
−
3
)
2
+
2
(
m
−
3
)
4
κ
2
+
(
m
−
3
)
2
−
4
κ
2
4
κ
(
m
−
3
)
−
m
−
3
κ
=
2
(
m
−
3
)
2
+
2
(
m
−
3
)
4
κ
2
+
(
m
−
3
)
2
4
κ
(
m
−
3
)
−
m
−
3
κ
=
(
m
−
3
)
+
4
κ
2
+
(
m
−
3
)
2
2
κ
−
2
(
m
−
3
)
2
κ
=
−
(
m
−
3
)
+
4
κ
2
+
(
m
−
3
)
2
2
κ
∈
(
0
,
1
)
1
−
b
0
1
+
b
0
=
m
−
3
+
2
κ
−
4
κ
2
+
(
m
−
3
)
2
m
−
3
−
2
κ
+
4
κ
2
+
(
m
−
3
)
2
=
[
m
−
3
+
2
κ
−
4
κ
2
+
(
m
−
3
)
2
]
[
(
m
−
3
−
2
κ
)
−
4
κ
2
+
(
m
−
3
)
2
]
[
(
m
−
3
−
2
κ
)
+
4
κ
2
+
(
m
−
3
)
2
]
[
(
m
−
3
−
2
κ
)
−
4
κ
2
+
(
m
−
3
)
2
]
=
[
(
m
−
3
−
4
κ
2
+
(
m
−
3
)
2
)
+
2
κ
]
[
(
m
−
3
−
4
κ
2
+
(
m
−
3
)
2
)
−
2
κ
]
(
m
−
3
−
2
κ
)
2
−
[
4
κ
2
+
(
m
−
3
)
2
]
=
(
m
−
3
−
4
κ
2
+
(
m
−
3
)
2
)
2
−
4
κ
2
[
(
m
−
3
)
2
+
4
κ
2
]
−
4
κ
(
m
−
3
)
−
[
4
κ
2
+
(
m
−
3
)
2
]
=
(
m
−
3
)
2
+
4
κ
2
+
(
m
−
3
)
2
−
2
(
m
−
3
)
4
κ
2
+
(
m
−
3
)
2
−
4
κ
2
−
4
κ
(
m
−
3
)
=
2
(
m
−
3
)
2
−
2
(
m
−
3
)
4
κ
2
+
(
m
−
3
)
2
−
4
κ
(
m
−
3
)
=
−
(
m
−
3
)
+
4
κ
2
+
(
m
−
3
)
2
2
κ
=
x
0
\begin{aligned} 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{aligned}
1+b01−b0x01+b01−b0========1+m−3−2κ+4κ2+(m−3)2=m−3m−3−2κ+4κ2+(m−3)2=m−3m−3+4κ2+(m−3)2−2κ=1−m−3−2κ+4κ2+(m−3)2=m−3m−3+2κ−4κ2+(m−3)2=m−3m−3+2κ−4κ2+(m−3)2=1−b01+b0−κm−3=m−3+2κ−4κ2+(m−3)2m−3−2κ+4κ2+(m−3)2−κm−3=(m−3+2κ)2−(4κ2+(m−3)2)(m−3+4κ2+(m−3)2−2κ)(m−3+4κ2+(m−3)2+2κ)−κm−3=(m−3)2+4κ2+4κ(m−3)−(4κ2+(m−3)2)(m−3+4κ2+(m−3)2)2−4κ2−κm−3=4κ(m−3)(m−3)2+4κ2+(m−3)2+2(m−3)4κ2+(m−3)2−4κ2−κm−3=4κ(m−3)2(m−3)2+2(m−3)4κ2+(m−3)2−κm−3=2κ(m−3)+4κ2+(m−3)2−2κ2(m−3)=2κ−(m−3)+4κ2+(m−3)2∈(0,1)m−3−2κ+4κ2+(m−3)2m−3+2κ−4κ2+(m−3)2[(m−3−2κ)+4κ2+(m−3)2][(m−3−2κ)−4κ2+(m−3)2][m−3+2κ−4κ2+(m−3)2][(m−3−2κ)−4κ2+(m−3)2](m−3−2κ)2−[4κ2+(m−3)2][(m−3−4κ2+(m−3)2)+2κ][(m−3−4κ2+(m−3)2)−2κ][(m−3)2+4κ2]−4κ(m−3)−[4κ2+(m−3)2](m−3−4κ2+(m−3)2)2−4κ2−4κ(m−3)(m−3)2+4κ2+(m−3)2−2(m−3)4κ2+(m−3)2−4κ2−4κ(m−3)2(m−3)2−2(m−3)4κ2+(m−3)22κ−(m−3)+4κ2+(m−3)2=x0 采样一个
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
=
1
−
(
1
+
b
0
)
y
1
−
(
1
−
b
0
)
y
∼
e
(
x
,
b
0
)
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)
y∼Be(2m−1,2m−1)⇒x=1−(1−b0)y1−(1+b0)y∼e(x,b0) 当
u
<
f
r
a
d
i
a
l
(
x
;
κ
,
m
)
M
∗
e
(
x
,
b
0
)
u \lt \frac{f_{radial}(x;\kappa,m)}{M*e(x,b_0)}
u<M∗e(x,b0)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
−
x
(
1
+
b
0
)
−
(
1
−
b
0
)
x
;
m
−
1
2
,
m
−
1
2
)
=
(
κ
/
2
)
ν
Γ
(
1
2
)
Γ
(
ν
+
1
2
)
I
ν
(
κ
)
B
(
m
−
1
2
,
m
−
1
2
)
b
0
m
−
3
2
[
(
1
+
b
0
)
−
(
1
−
b
0
)
x
]
m
−
3
e
x
p
(
κ
x
)
/
M
=
[
(
1
+
b
0
)
−
(
1
−
b
0
)
x
]
m
−
3
e
x
p
(
κ
x
)
[
(
1
+
b
0
)
−
(
1
−
b
0
)
x
0
]
m
−
3
e
x
p
(
κ
x
0
)
=
(
(
1
+
b
0
)
−
(
1
−
b
0
)
x
(
1
+
b
0
)
−
(
1
−
b
0
)
x
0
)
m
−
3
e
x
p
(
κ
(
x
−
x
0
)
)
l
o
g
f
r
a
d
i
a
l
(
x
;
κ
,
m
)
M
∗
e
(
x
,
b
0
)
=
(
m
−
3
)
l
o
g
(
(
1
+
b
0
)
−
(
1
−
b
0
)
x
(
1
+
b
0
)
−
(
1
−
b
0
)
x
0
)
+
κ
(
x
−
x
0
)
=
(
m
−
3
)
l
o
g
(
1
−
1
−
b
0
1
+
b
0
x
1
−
1
−
b
0
1
+
b
0
x
0
)
+
κ
(
x
−
x
0
)
=
(
m
−
3
)
l
o
g
(
1
−
x
0
x
1
−
x
0
2
)
+
κ
(
x
−
x
0
)
=
(
m
−
3
)
l
o
g
(
1
−
x
0
x
)
−
(
m
−
3
)
l
o
g
(
1
−
x
0
2
)
+
κ
(
x
−
x
0
)
=
κ
x
+
(
m
−
3
)
l
o
g
(
1
−
x
0
x
)
−
[
κ
x
0
+
(
m
−
3
)
l
o
g
(
1
−
x
0
2
)
]
=
κ
x
+
(
m
−
3
)
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-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{\Beta(\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*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_0x}{1-x_0^2}\right) + \kappa (x - x_0) \\ &= (m-3)log(1-x_0x) - (m-3)log(1-x_0^2) + \kappa (x - x_0) \\ &= \kappa x + (m-3)log(1-x_0x) - [\kappa x_0 + (m-3)log(1-x_0^2)] \\ &= \kappa x + (m-3)log(1-x_0x) - c \end{aligned}
M∗e(x,b0)fradial(x;κ,m)logM∗e(x,b0)fradial(x;κ,m)=Mf((1+b0)−(1−b0)x1−x;2m−1,2m−1)fradial(x;κ,m)=Γ(21)Γ(ν+21)Iν(κ)(κ/2)νb02m−3B(2m−1,2m−1)[(1+b0)−(1−b0)x]m−3exp(κx)/M=[(1+b0)−(1−b0)x0]m−3exp(κx0)[(1+b0)−(1−b0)x]m−3exp(κx)=((1+b0)−(1−b0)x0(1+b0)−(1−b0)x)m−3exp(κ(x−x0))=(m−3)log((1+b0)−(1−b0)x0(1+b0)−(1−b0)x)+κ(x−x0)=(m−3)log(1−1+b01−b0x01−1+b01−b0x)+κ(x−x0)=(m−3)log(1−x021−x0x)+κ(x−x0)=(m−3)log(1−x0x)−(m−3)log(1−x02)+κ(x−x0)=κx+(m−3)log(1−x0x)−[κx0+(m−3)log(1−x02)]=κx+(m−3)log(1−x0x)−c 即, 当
l
o
g
u
<
κ
x
+
(
m
−
3
)
l
o
g
(
1
−
x
0
x
)
−
c
logu < \kappa x + (m-3)log(1-x_0x) - c
logu<κx+(m−3)log(1−x0x)−c 时接受样本
x
x
x. 与 Communications in Statistics - Simulation and Computation [2] 中的判别公式一致了(除了
(
m
−
3
)
(m-3)
(m−3) 和
(
m
−
1
)
(m-1)
(m−1)):
但是与 [1] 不一致, 且 [2] 声称 [1] 的采样方案有问题:
意思大概是作者在实践中并没有的到正确的采样, 且不知道问题出在哪里, 所以他们又根据 [1] 推导了一遍, 得出了上述接受-拒绝判别式. 但, 我也推导了一遍, 发现两者是等价的…
问题到底在哪里?
似乎所有的差别就在
e
(
x
,
b
)
e(x, b)
e(x,b), 到底是:
e
(
x
,
b
)
=
b
m
−
3
2
B
(
m
−
1
2
,
m
−
1
2
)
(
1
−
x
2
)
m
−
3
2
[
(
1
+
b
)
−
(
1
−
b
)
x
]
m
−
3
\begin{aligned} e(x,b) =& \frac{b^{\frac{m-3}{2}}}{\Beta(\frac{m-1}{2},\frac{m-1}{2})} \frac{(1-x^2)^{\frac{m-3}{2}}}{[(1+b)-(1-b)x]^{m-3}} \end{aligned}
e(x,b)=B(2m−1,2m−1)b2m−3[(1+b)−(1−b)x]m−3(1−x2)2m−3 还是:
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
\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}} \end{aligned}
e(x,b)=B(2m−1,2m−1)2b2m−1[(1+b)−(1−b)x]m−1(1−x2)2m−3 ? ? ? ? ? ? ? 找到问题出在哪了, 保留以上过程, 以警钟长鸣 ? ? ? ? ? ? ? ?
记住, 永远不要想当然! 概率密度函数的变换不是进行一下变量代换就 OK 了的! 前面
t
t
t 的概率密度函数推导中, 就知道
t
t
t 和
θ
\theta
θ 在变换的时候需要带上
d
t
,
d
θ
dt, d\theta
dt,dθ, 怎么现在就忘了? 根本原因是对知识的理解不够透彻. 我们的研究对象本质上是概率
f
(
y
;
α
,
β
)
d
y
f(y; \alpha, \beta)dy
f(y;α,β)dy, 而不是概率密度
f
(
y
;
α
,
β
)
f(y; \alpha, \beta)
f(y;α,β), 当进行变量代换时, 不光要变换概率密度中的变量,
d
y
dy
dy 也要变换, 甚至求
[
a
,
b
]
[a, b]
[a,b] 区间上的概率时,
∫
a
b
f
(
y
;
α
,
β
)
d
y
\int_a^b f(y; \alpha, \beta) dy
∫abf(y;α,β)dy, 积分区间也要变换.
在得到随机变量变换关系: X = 1 − ( 1 + b ) Y 1 − ( 1 − b ) Y ⇔ Y = 1 − X ( 1 + b ) − ( 1 − b ) X \begin{aligned} X = \frac{1-(1+b)Y}{1-(1-b)Y} &\Leftrightarrow Y = \frac{1-X}{(1+b)-(1-b)X} \end{aligned} X=1−(1−b)Y1−(1+b)Y⇔Y=(1+b)−(1−b)X1−X 后, 假如让你积分求区间 [ y 1 , y 2 ] [y_1,y_2] [y1,y2] 上的概率呢? ∫ y 1 y 2 f ( y ; α , β ) d y = ∫ 1 − x 1 ( 1 + b ) − ( 1 − b ) x 1 1 − x 2 ( 1 + b ) − ( 1 − b ) x 2 f ( 1 − x ( 1 + b ) − ( 1 − b ) x ; α , β ) d ( 1 − x ( 1 + b ) − ( 1 − b ) x ) = ∫ 1 − x 1 ( 1 + b ) − ( 1 − b ) x 1 1 − x 2 ( 1 + b ) − ( 1 − b ) x 2 f ( 1 − x ( 1 + b ) − ( 1 − b ) x ; α , β ) − 2 b [ ( 1 + b ) − ( 1 − b ) x ] 2 d x \begin{aligned} \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{aligned} ∫y1y2f(y;α,β)dy=∫(1+b)−(1−b)x11−x1(1+b)−(1−b)x21−x2f((1+b)−(1−b)x1−x;α,β)d((1+b)−(1−b)x1−x)=∫(1+b)−(1−b)x11−x1(1+b)−(1−b)x21−x2f((1+b)−(1−b)x1−x;α,β)[(1+b)−(1−b)x]2−2bdx
d y = d ( 1 − x ( 1 + b ) − ( 1 − b ) x ) = − [ ( 1 + b ) − ( 1 − b ) x ] − ( 1 − x ) [ − ( 1 − b ) ] [ ( 1 + b ) − ( 1 − b ) x ] 2 d x = − ( 1 + b ) + ( 1 − b ) x + ( 1 − b ) − ( 1 − b ) x [ ( 1 + b ) − ( 1 − b ) x ] 2 d x = − 2 b [ ( 1 + b ) − ( 1 − b ) x ] 2 d x \begin{aligned} 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{aligned} dy=d((1+b)−(1−b)x1−x)=[(1+b)−(1−b)x]2−[(1+b)−(1−b)x]−(1−x)[−(1−b)]dx=[(1+b)−(1−b)x]2−(1+b)+(1−b)x+(1−b)−(1−b)xdx=[(1+b)−(1−b)x]2−2bdx
由于 − 2 b [ ( 1 + b ) − ( 1 − b ) x ] 2 < 0 \frac{-2b}{[(1+b)-(1-b)x]^2} \lt 0 [(1+b)−(1−b)x]2−2b<0, 变量 x , y x, y x,y 之间是负相关的, 那么积分式写为: ∫ 1 − x 2 ( 1 + b ) − ( 1 − b ) x 2 1 − x 1 ( 1 + b ) − ( 1 − b ) x 1 f ( 1 − x ( 1 + b ) − ( 1 − b ) x ; α , β ) 2 b [ ( 1 + b ) − ( 1 − b ) x ] 2 d x \begin{aligned} \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{aligned} ∫(1+b)−(1−b)x21−x2(1+b)−(1−b)x11−x1f((1+b)−(1−b)x1−x;α,β)[(1+b)−(1−b)x]22bdx 所以, 关于 x x x 的概率密度函数是: e ( x , b ) = f ( 1 − x ( 1 + b ) − ( 1 − b ) x ; α , β ) 2 b [ ( 1 + b ) − ( 1 − b ) x ] 2 = b m − 3 2 B ( m − 1 2 , m − 1 2 ) ( 1 − x 2 ) m − 3 2 [ ( 1 + b ) − ( 1 − b ) x ] m − 3 2 b [ ( 1 + b ) − ( 1 − b ) x ] 2 = 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 \begin{aligned} 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}}}{\Beta(\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}}}{\Beta(\frac{m-1}{2},\frac{m-1}{2})} \frac{(1-x^2)^{\frac{m-3}{2}}}{[(1+b)-(1-b)x]^{m-1}} \end{aligned} e(x,b)===f((1+b)−(1−b)x1−x;α,β)[(1+b)−(1−b)x]22bB(2m−1,2m−1)b2m−3[(1+b)−(1−b)x]m−3(1−x2)2m−3[(1+b)−(1−b)x]22bB(2m−1,2m−1)2b2m−1[(1+b)−(1−b)x]m−1(1−x2)2m−3 这下一致了.
平均接受率
在拒绝采样法的简介中说: 平均接受率为
1
c
\frac{1}{c}
c1, 那么本拒绝采样法的平均接受率为
1
M
\frac{1}{M}
M1.
1
M
=
e
(
x
0
,
b
0
)
f
(
x
0
,
κ
,
m
)
\begin{aligned} \frac{1}{M} = \frac{e(x_0, b_0)}{f(x_0, \kappa, m)} \end{aligned}
M1=f(x0,κ,m)e(x0,b0) 公式太过复杂, 简化无果后, 在 desmos 网站画出了其随
κ
,
m
\kappa, m
κ,m 变化的图像(计算比较复杂,可能需要等一会儿图像才会出来):
可见, 接受率还是比较高的, 趋势大概是:
m
m
m 越大, 接受率越高,
κ
\kappa
κ 越大, 接受率越低. 但即使在
m
=
4
,
κ
=
20
m=4, \kappa=20
m=4,κ=20 时, 接受率依然大于
0.7
0.7
0.7.
下面给出了
(
m
=
50
,
κ
=
5
)
(m=50,\kappa=5)
(m=50,κ=5) 和
(
m
=
5
,
κ
=
20
)
(m=5,\kappa=20)
(m=5,κ=20) 时,
f
r
a
d
i
a
l
(
x
;
κ
,
m
)
f_{radial}(x;\kappa,m)
fradial(x;κ,m) 和
M
∗
e
(
x
,
b
0
)
M*e(x, b_0)
M∗e(x,b0) 的图像(前者为蓝色, 后者为黑色). 可见, 当
(
m
=
50
,
κ
=
5
)
(m=50,\kappa=5)
(m=50,κ=5) 时, 两个图像几乎是重合的; 当
(
m
=
5
,
κ
=
20
)
(m=5,\kappa=20)
(m=5,κ=20) 时, 两者出现了较大出入, 表现为接受率的下降.
参考文献
[1] Computer Generation of Distributions on the m-spher.
[2] Communications in Statistics - Simulation and Computation.