PGM模型的典型问题之一就是推断问题。
推断
推断问题定义为，已知模型中的部分随机变量 e \textbf e e，计算其他部分随机变量 q \textbf q q的后验概率，即:
P ( q ∣ e ) P(\textbf q|\textbf e) P(q∣e)
一般的根据贝叶斯公式直接进行计算:
P ( q ∣ e ) = P ( q , e ) P ( e ) P(\textbf q|\textbf e)=\frac{P(\textbf q,\textbf e)}{P(\textbf e)} P(q∣e)=P(e)P(q,e)
但是在模型比较复杂，涉及的变量数量庞大的时候，直接精确计算开销太大，一般会采用近似计算的方法。采样法是比较常用的一种近似计算的方法，也叫蒙特卡罗法(Monte Carlo Method)。
蒙特卡罗法
采样法是采集符合 p ( x ) p(x) p(x)分布的一定数量的样本，并通过样本来估计此分布的一些值，比如期望值等。
例如，随机变量 x x x的分布函数为 p ( x ) ， 为 了 计 算 函 数 f ( x ) 的 期 望 p(x)，为了计算函数f(x)的期望 p(x)，为了计算函数f(x)的期望：
E [ f ( x ) ] = ∫ f ( x ) p ( x ) d x E[f(x)]=\int_{}f(x)p(x)dx E[f(x)]=∫f(x)p(x)dx
当 f ( x ) p ( x ) f(x)p(x) f(x)p(x)的形式比较复杂，不适合直接数值计算的情况下，可有采用采样法，计算一个估计值 f ^ ( x ) \hat f(x) f^(x)：
f ^ ( x ) = 1 N ( f ( x ( 1 ) ) + f ( x ( 2 ) ) + . . . + f ( x ( N ) ) ) \hat f(x)=\frac{1}{N}(f(x^{(1)})+f(x^{(2)})+...+f(x^{(N)})) f^(x)=N1(f(x(1))+f(x(2))+...+f(x(N)))
f ^ ( x ) → E [ f ( x ) ] , N → ∞ \hat f(x)\rightarrow E[f(x)],N\rightarrow\infty f^(x)→E[f(x)],N→∞
其中 x ( 1 ) , x ( 2 ) , . . . , x ( N ) x^{(1)},x^{(2)},...,x^{(N)} x(1),x(2),...,x(N)为符合分布函数 p ( x ) p(x) p(x)的一批样本。
采样法的难点在于，如何生成符合概率分布p(x)的样本。 x ( 1 ) , x ( 2 ) , . . . , x ( N ) x^{(1)},x^{(2)},...,x^{(N)} x(1),x(2),...,x(N)。
我们知道，计算机比较容易的生成均匀分布的样本，如何生成符合给定概率密度函数 p ( x ) p(x) p(x)的样本呢？
假定随机变量 x x x符合概率密度函数 p ( x ) p(x) p(x)， x x x的累积概率函数 c ( x ) c(x) c(x)；另外假设随机变量 γ \gamma γ在 [ 0 , 1 ] [0,1] [0,1]区间内均匀分布，即:
p
d
f
(
x
)
=
p
(
x
)
pdf(x)=p(x)
pdf(x)=p(x)
c
d
f
(
x
)
=
c
(
x
)
cdf(x)=c(x)
cdf(x)=c(x)
p
d
f
(
γ
)
=
1
pdf(\gamma)=1
pdf(γ)=1
c
d
f
(
γ
)
=
γ
cdf(\gamma)=\gamma
cdf(γ)=γ
将 c ( x ) c(x) c(x)的逆函数记为 c − 1 ( x ) c^{-1}(x) c−1(x). 则随机变量 x = c − 1 ( γ ) x=c^{-1}(\gamma) x=c−1(γ)的累积概率函数为 c ( x ) c(x) c(x),概率密度函数为 p ( x ) p(x) p(x).
由累积概率函数的定义开始，简单证明如下：
c
d
f
(
x
)
=
P
(
c
−
1
(
γ
)
≤
x
)
cdf(x)=P(c^{-1}(\gamma)\le x)
cdf(x)=P(c−1(γ)≤x)
=
P
(
γ
≤
c
(
x
)
)
=P(\gamma\le c(x))
=P(γ≤c(x))
=
c
(
x
)
=c(x)
=c(x)
当 p ( x ) p(x) p(x)的形式比较复杂， c ( x ) c(x) c(x)的逆函数很难计算的情况下，可以采用拒绝采样，重要性采样，马尔科夫蒙特卡洛采样。
拒绝采样
如果分布 p ( x ) p(x) p(x)不能直接采样，可以对一个提议分布 q ( x ) q(x) q(x)进行采用，然后按照一定的标准拒绝一部分样本，使最终的样本符合分布 p ( x ) p(x) p(x).
例如对于分布 p ( x ) p(x) p(x),提议分布 q ( x ) q(x) q(x)和常数 k k k满足 k q ( x ) ≥ p ( x ) kq(x)\ge p(x) kq(x)≥p(x),并以以下接受概率接受样本：
α ( x ) = p ( x ) k q ( x ) \alpha(x)=\frac{p(x)}{kq(x)} α(x)=kq(x)p(x)
重要性采样
E
p
[
f
(
x
)
]
=
∫
x
f
(
x
)
p
(
x
)
d
x
E_p[f(x)]=\int_xf(x)p(x)dx
Ep[f(x)]=∫xf(x)p(x)dx
=
∫
x
f
(
x
)
p
(
x
)
q
(
x
)
q
(
x
)
d
x
=\int_xf(x)\frac{p(x)}{q(x)}q(x)dx
=∫xf(x)q(x)p(x)q(x)dx
=
∫
x
f
(
x
)
w
(
x
)
q
(
x
)
d
x
=\int_xf(x)w(x)q(x)dx
=∫xf(x)w(x)q(x)dx
=
E
q
[
f
(
x
)
w
(
x
)
]
=E_q[f(x)w(x)]
=Eq[f(x)w(x)]
所以：
E
p
(
x
)
=
E
q
(
f
(
x
)
w
(
x
)
)
E_p(x)=E_q(f(x)w(x))
Ep(x)=Eq(f(x)w(x))
≈
1
N
(
f
(
x
(
1
)
)
w
(
x
(
1
)
)
+
f
(
x
(
2
)
)
w
(
x
(
2
)
)
+
.
.
.
+
f
(
x
(
N
)
)
w
(
x
(
N
)
)
\approx \frac{1}{N}(f(x^{(1)})w(x^{(1)})+f(x^{(2)})w(x^{(2)})+...+f(x^{(N)})w(x^{(N)})
≈N1(f(x(1))w(x(1))+f(x(2))w(x(2))+...+f(x(N))w(x(N))
其中 x ( 1 ) , x ( 2 ) , . . . , x ( N ) x^{(1)},x^{(2)},...,x^{(N)} x(1),x(2),...,x(N)为符合分布q(x)的样本， w ( x ( i ) ) w(x^{(i)}) w(x(i))样本的权重。
马尔科夫蒙特卡洛方法
MCMC方法的关键是构造一个平稳分布为p(x)的马尔可夫链，并且马尔可夫链的转移概率 q ( x t ∣ x t − 1 ) q(x_t|x_{t-1}) q(xt∣xt−1)比较容易采用，然后从达到平稳状态的马尔可夫链中抽取样本 x x x，这样样本 x x x就符合分布 p ( x ) p(x) p(x).
MCMC方法有很多中实现算法，比如MH算法为例：MH采用的马尔科夫状态转移分布为 q ( x t ∣ x t − 1 ) q(x_t|x_{t-1}) q(xt∣xt−1)，这样的马尔科夫链平稳状态的分布往往不会是 p ( x ) p(x) p(x),MH算法中采用拒绝采样的思想，以一定的概率接受样本,接受概率为
α ( x t , x t − 1 ) = m i n ( 1 , p ( x t ) q ( x t − 1 ∣ x t ) p ( x t − 1 ) q ( x t ∣ x t − 1 ) ) \alpha(x_t,x_{t-1})=min(1,\frac{p(x_t)q(x_{t-1}|x_t)}{p(x_{t-1})q(x_t|x_{t-1})}) α(xt,xt−1)=min(1,p(xt−1)q(xt∣xt−1)p(xt)q(xt−1∣xt))
所以实际的状态转移分布为:
q ′ ( x t ∣ x t − 1 ) = q ( x t ∣ x t − 1 ) α ( x t , x t − 1 ) q'(x_t|x_{t-1})=q(x_t|x_{t-1})\alpha(x_t,x_{t-1}) q′(xt∣xt−1)=q(xt∣xt−1)α(xt,xt−1)
状态转移分布为 q ′ ( x t ∣ x t − 1 ) q'(x_t|x_{t-1}) q′(xt∣xt−1)的马尔可夫链的平稳分布为p(x).
简单证明如下：
p
(
x
t
−
1
)
q
′
(
x
t
∣
x
t
−
1
)
=
p
(
x
t
−
1
)
q
(
x
t
∣
x
t
−
1
)
α
(
x
t
,
x
t
−
1
)
p(x_{t-1})q'(x_t|x_{t-1})=p(x_{t-1})q(x_t|x_{t-1})\alpha(x_t,x_{t-1})
p(xt−1)q′(xt∣xt−1)=p(xt−1)q(xt∣xt−1)α(xt,xt−1)
=
p
(
x
t
−
1
)
q
(
x
t
∣
x
t
−
1
)
m
i
n
(
1
,
p
(
x
t
)
q
(
x
t
−
1
∣
x
t
)
p
(
x
t
−
1
)
q
(
x
t
∣
x
t
−
1
)
)
)
=p(x_{t-1})q(x_t|x_{t-1})min(1,\frac{p(x_t)q(x_{t-1}|x_t)}{p(x_{t-1})q(x_t|x_{t-1})}))
=p(xt−1)q(xt∣xt−1)min(1,p(xt−1)q(xt∣xt−1)p(xt)q(xt−1∣xt)))
=
m
i
n
(
p
(
x
t
−
1
)
q
(
x
t
∣
x
t
−
1
)
,
p
(
x
t
)
q
(
x
t
−
1
∣
x
t
)
)
=min(p(x_{t-1})q(x_t|x_{t-1}),p(x_t)q(x_{t-1}|x_t))
=min(p(xt−1)q(xt∣xt−1),p(xt)q(xt−1∣xt))
=
p
(
x
t
)
q
(
x
t
−
1
∣
x
t
)
m
i
n
(
p
(
x
t
−
1
)
q
(
x
t
∣
x
t
−
1
)
p
(
x
t
)
q
(
x
t
−
1
∣
x
t
)
,
1
)
=p(x_t)q(x_{t-1}|x_t)min(\frac{p(x_{t-1})q(x_t|x_{t-1})}{p(x_t)q(x_{t-1}|x_t)}, 1)
=p(xt)q(xt−1∣xt)min(p(xt)q(xt−1∣xt)p(xt−1)q(xt∣xt−1),1)
=
p
(
x
t
)
q
(
x
t
−
1
∣
x
t
)
α
(
x
t
−
1
,
x
t
)
=p(x_t)q(x_{t-1}|x_t)\alpha(x_{t-1},x_t)
=p(xt)q(xt−1∣xt)α(xt−1,xt)
=
p
(
x
t
)
q
′
(
x
t
−
1
∣
x
t
)
=p(x_t)q'(x_{t-1}|x_t)
=p(xt)q′(xt−1∣xt)
根据马尔科夫链的细致平衡条件，此马尔可夫链的稳定状态的分布为 p ( x ) p(x) p(x)
扫一扫
7632
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
