在前面《MCMC(一):蒙特卡罗方法和马尔科夫链》中,我们讲到给定一个概率平稳分布 π \pi π, 很难直接找到对应的马尔科夫链状态转移矩阵 P P P。而只要解决这个问题,我们就可以找到一种通用的概率分布采样方法,进而用于蒙特卡罗模拟。本篇我们就讨论解决这个问题的办法:MCMC采样和它的易用版M-H采样。
本篇博客主要转自参考文献【1】,在原文的基础上,为了更容易增加理解,略有删改增。
一、马尔科夫链的细致平稳条件
在解决从平稳分布 π \pi π, 找到对应的马尔科夫链状态转移矩阵 P P P之前,我们还需要先看看马尔科夫链的细致平稳条件。定义如下:
如果非周期马尔科夫链的状态转移矩阵 P P P和概率分布 π ( x ) \pi(x) π(x)对于所有的 i , j i,j i,j满足: π ( i ) P ( i , j ) = π ( j ) P ( j , i ) \pi(i)P(i,j) = \pi(j)P(j,i) π(i)P(i,j)=π(j)P(j,i)
则称概率分布 π ( x ) \pi(x) π(x)是状态转移矩阵 P P P的平稳分布。
证明很简单,由细致平稳条件有:
∑
i
=
1
∞
π
(
i
)
P
(
i
,
j
)
=
∑
i
=
1
∞
π
(
j
)
P
(
j
,
i
)
=
π
(
j
)
∑
i
=
1
∞
P
(
j
,
i
)
=
π
(
j
)
\sum\limits_{i=1}^{\infty}\pi(i)P(i,j) = \sum\limits_{i=1}^{\infty} \pi(j)P(j,i) = \pi(j)\sum\limits_{i=1}^{\infty} P(j,i) = \pi(j)
i=1∑∞π(i)P(i,j)=i=1∑∞π(j)P(j,i)=π(j)i=1∑∞P(j,i)=π(j)
将上式用矩阵表示即为:
π
P
=
π
\pi P = \pi
πP=π
即满足马尔可夫链的收敛性质。也就是说,只要我们找到了可以使概率分布 π ( x ) \pi(x) π(x)满足细致平稳分布的矩阵 P P P即可。这给了我们寻找从平稳分布 π \pi π, 找到对应的马尔科夫链状态转移矩阵 P P P的新思路。
不过不幸的是,仅仅从细致平稳条件还是很难找到合适的矩阵
P
P
P。比如我们的目标平稳分布是
π
(
x
)
\pi(x)
π(x),随机找一个马尔科夫链状态转移矩阵
Q
Q
Q,它是很难满足细致平稳条件的,即:
π
(
i
)
Q
(
i
,
j
)
≠
π
(
j
)
Q
(
j
,
i
)
\pi(i)Q(i,j) \neq \pi(j)Q(j,i)
π(i)Q(i,j)̸=π(j)Q(j,i)
那么如何使这个等式满足呢?下面我们来看MCMC采样如何解决这个问题。
二、MCMC采样
由于一般情况下,目标平稳分布
π
(
x
)
\pi(x)
π(x)和某一个马尔科夫链状态转移矩阵
Q
Q
Q不满足细致平稳条件,即
π
(
i
)
Q
(
i
,
j
)
≠
π
(
j
)
Q
(
j
,
i
)
\pi(i)Q(i,j) \neq \pi(j)Q(j,i)
π(i)Q(i,j)̸=π(j)Q(j,i)
我们可以对上式做一个改造,使细致平稳条件成立。方法是引入一个
α
(
i
,
j
)
\alpha(i,j)
α(i,j),使上式可以取等号,即:
π
(
i
)
Q
(
i
,
j
)
α
(
i
,
j
)
=
π
(
j
)
Q
(
j
,
i
)
α
(
j
,
i
)
\pi(i)Q(i,j)\alpha(i,j) = \pi(j)Q(j,i)\alpha(j,i)
π(i)Q(i,j)α(i,j)=π(j)Q(j,i)α(j,i)
问题是什么样的
α
(
i
,
j
)
\alpha(i,j)
α(i,j)可以使等式成立呢?其实很简单,只要满足下两式即可:
α
(
i
,
j
)
=
π
(
j
)
Q
(
j
,
i
)
\alpha(i,j) = \pi(j)Q(j,i)
α(i,j)=π(j)Q(j,i)
α ( j , i ) = π ( i ) Q ( i , j ) \alpha(j,i) = \pi(i)Q(i,j) α(j,i)=π(i)Q(i,j)
这样,我们就得到了我们的分布
π
(
x
)
\pi(x)
π(x)对应的马尔科夫链状态转移矩阵
P
P
P,满足:
P
(
i
,
j
)
=
Q
(
i
,
j
)
α
(
i
,
j
)
P(i,j) = Q(i,j)\alpha(i,j)
P(i,j)=Q(i,j)α(i,j)
也就是说,我们的目标矩阵 P P P可以通过任意一个马尔科夫链状态转移矩阵 Q Q Q乘以 α ( i , j ) \alpha(i,j) α(i,j)得到。 α ( i , j ) \alpha(i,j) α(i,j)我们有一般称之为接受率,取值在 [ 0 , 1 ] [0,1] [0,1]之间。物理意义可以理解为:在原来的马氏链上,从状态 i i i以 Q ( i , j ) Q(i,j) Q(i,j)的概率转移到状态 j j j的时候,我们以 α ( i , j ) \alpha(i,j) α(i,j)的概率接受这个转移。 即目标矩阵 P P P可以通过任意一个马尔科夫链状态转移矩阵 Q Q Q以一定的接受率获得。这个很像拒绝采样的思想(详见《拒绝采样(Reject Sampling)原理详解》),那里是以一个常用分布通过一定的接受-拒绝概率得到一个非常见分布,这里是以一个常见的马尔科夫链状态转移矩阵 Q Q Q通过一定的接受-拒绝概率得到目标转移矩阵 P P P,两者的解决问题思路是类似的。
好了,现在我们来总结下MCMC的采样过程。
1)输入我们任意选定的马尔科夫链状态转移矩阵 Q Q Q,平稳分布 π ( x ) \pi(x) π(x),设定状态转移次数阈值 n 1 n_1 n1,需要的样本个数 n 2 n_2 n2
2)从任意简单概率分布采样得到初始状态值 x 0 x_0 x0
3)for t = 0 t=0 t=0 to n 1 + n 2 − 1 n_1+n_2−1 n1+n2−1:
a) 从条件概率分布
Q
(
x
∣
x
t
)
Q(x|x_t)
Q(x∣xt)中采样得到样本
x
∗
x_∗
x∗
b) 从均匀分布采样
u
∼
u
n
i
f
o
r
m
[
0
,
1
]
u \sim uniform[0,1]
u∼uniform[0,1]
c) 如果
u
<
α
(
x
t
,
x
∗
)
=
π
(
x
∗
)
Q
(
x
∗
,
x
t
)
u < \alpha(x_t,x_{*}) = \pi(x_{*})Q(x_{*},x_t)
u<α(xt,x∗)=π(x∗)Q(x∗,xt), 则接受转移
x
t
→
x
∗
x_t \to x_{*}
xt→x∗,即
x
t
+
1
=
x
∗
x_{t+1}= x_{*}
xt+1=x∗
d) 否则不接受转移,即
x
t
+
1
=
x
t
x_{t+1}= x_{t}
xt+1=xt
样本集 ( x n 1 , x n 1 + 1 , . . . , x n 1 + n 2 − 1 ) (x_{n_1}, x_{n_1+1},..., x_{n_1+n_2-1}) (xn1,xn1+1,...,xn1+n2−1)即为我们需要的平稳分布对应的样本集。对4-d多说一点: x t + 1 = x t x_{t+1}= x_{t} xt+1=xt这意味着,尽管不接受转移,但是也把它放入了样本序列。而在有些算法中,被拒绝的样本就不进入采样样本序列。尽管双方有些mismatch,不过对理解算法影响不大。
上面这个过程基本上就是MCMC采样的完整采样理论了,但是这个采样算法还是比较难在实际中应用,为什么呢?问题在上面第三步的c步骤,接受率这儿。由于 α ( x t , x ∗ ) \alpha(x_t,x_{*}) α(xt,x∗)可能非常的小,比如0.1,导致我们大部分的采样值都被拒绝转移,采样效率很低。有可能我们采样了上百万次马尔可夫链还没有收敛,也就是上面这个 n 1 n_1 n1要非常非常的大,这让人难以接受,怎么办呢?这时就轮到我们的M-H采样出场了。
三、M-H采样
M-H采样是Metropolis-Hastings采样的简称,这个算法首先由Metropolis提出,被Hastings改进,因此被称之为Metropolis-Hastings采样或M-H采样
M-H采样解决了我们上一节MCMC采样接受率过低的问题。
我们回到MCMC采样的细致平稳条件:
π
(
i
)
Q
(
i
,
j
)
α
(
i
,
j
)
=
π
(
j
)
Q
(
j
,
i
)
α
(
j
,
i
)
\pi(i)Q(i,j)\alpha(i,j) = \pi(j)Q(j,i)\alpha(j,i)
π(i)Q(i,j)α(i,j)=π(j)Q(j,i)α(j,i)
我们采样效率低的原因是
α
(
i
,
j
)
\alpha(i,j)
α(i,j)太小了,比如为0.1,而
α
(
i
,
j
)
\alpha(i,j)
α(i,j)为0.2。即:
π
(
i
)
Q
(
i
,
j
)
×
0.1
=
π
(
j
)
Q
(
j
,
i
)
×
0.2
\pi(i)Q(i,j)\times 0.1 = \pi(j)Q(j,i)\times 0.2
π(i)Q(i,j)×0.1=π(j)Q(j,i)×0.2
这时我们可以看到,如果两边同时扩大五倍,接受率提高到了0.5,但是细致平稳条件却仍然是满足的,即:
π
(
i
)
Q
(
i
,
j
)
×
0.5
=
π
(
j
)
Q
(
j
,
i
)
×
1
\pi(i)Q(i,j)\times 0.5 = \pi(j)Q(j,i)\times 1
π(i)Q(i,j)×0.5=π(j)Q(j,i)×1
这样我们的接受率可以做如下改进,即:
α
(
i
,
j
)
=
m
i
n
{
π
(
j
)
Q
(
j
,
i
)
π
(
i
)
Q
(
i
,
j
)
,
1
}
\alpha(i,j) = min\{ \frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)},1\}
α(i,j)=min{π(i)Q(i,j)π(j)Q(j,i),1}
对上式多解释一下:此时 α ( i , j ) \alpha(i,j) α(i,j)已经不是原来的意义,变为了 α ( i , j ) / α ( j , i ) \alpha(i,j)/\alpha(j,i) α(i,j)/α(j,i)的倍率关系,这里默认右边为1,看左边的较大的那个接受率的值。即当原始的 α ( i , j ) < α ( j , i ) \alpha(i,j)<\alpha(j,i) α(i,j)<α(j,i)时, α ( i , j ) \alpha(i,j) α(i,j)扩大了 1 / α ( j , i ) 1/\alpha(j,i) 1/α(j,i)倍;当原始的 α ( i , j ) > α ( j , i ) \alpha(i,j)>\alpha(j,i) α(i,j)>α(j,i)时, α ( i , j ) \alpha(i,j) α(i,j)直接变为1。
通过这个微小的改造,我们就得到了可以在实际应用中使用的M-H采样算法过程如下:
1)输入我们任意选定的马尔科夫链状态转移矩阵 Q Q Q,平稳分布 π ( x ) \pi(x) π(x),设定状态转移次数阈值 n 1 n_1 n1,需要的样本个数 n 2 n_2 n2
2)从任意简单概率分布采样得到初始状态值 x 0 x_0 x0
3)for t = 0 t=0 t=0 to n 1 + n 2 − 1 n_1+n_2−1 n1+n2−1:
a) 从条件概率分布
Q
(
x
∣
x
t
)
Q(x|x_t)
Q(x∣xt)中采样得到样本
x
∗
x_∗
x∗
b) 从均匀分布采样
u
∼
u
n
i
f
o
r
m
[
0
,
1
]
u \sim uniform[0,1]
u∼uniform[0,1]
c) 如果
u
<
α
(
x
t
,
x
∗
)
=
m
i
n
{
π
(
j
)
Q
(
j
,
i
)
π
(
i
)
Q
(
i
,
j
)
,
1
}
u < \alpha(x_t,x_{*}) = min\{ \frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)},1\}
u<α(xt,x∗)=min{π(i)Q(i,j)π(j)Q(j,i),1}, 则接受转移
x
t
→
x
∗
x_t \to x_{*}
xt→x∗,即
x
t
+
1
=
x
∗
x_{t+1}= x_{*}
xt+1=x∗
d) 否则不接受转移,即
x
t
+
1
=
x
t
x_{t+1}= x_{t}
xt+1=xt
样本集 ( x n 1 , x n 1 + 1 , . . . , x n 1 + n 2 − 1 ) (x_{n_1}, x_{n_1+1},..., x_{n_1+n_2-1}) (xn1,xn1+1,...,xn1+n2−1)即为我们需要的平稳分布对应的样本集。
四、M-H采样实例
为了更容易理解,这里给出一个M-H采样的实例。
在例子里,我们的目标平稳分布是一个均值3,标准差2的正态分布,而选择的马尔可夫链状态转移矩阵 Q ( i , j ) Q(i,j) Q(i,j)的条件转移概率是以 i i i为均值,方差1的正态分布在位置 j j j的值。
注:这个例子仅仅用来让大家加深对M-H采样过程的理解。毕竟一个普通的一维正态分布用不着去用M-H采样来获得样本。
代码如下:
import random
from scipy.stats import norm
import matplotlib.pyplot as plt
def norm_dist_prob(theta):
#求概率密度函数指定点(theta)的函数值,即theta点附近的可能性。
y = norm.pdf(theta, loc=3, scale=2)
return y
# n1+n2=5000
T = 5000
# 初始状态值x_0就是pi[0],即为0。
pi = [0 for i in range(T)]
t = 0
while t < T-1:
# pi_star可以认为是第j个状态,pi[t - 1]可以认为是第i个状态。
t = t + 1
# 生成服从均值是pi[t - 1],标准差是1的正态分布的随机数。pi_star就是x_star。
# 也就是说pi[t - 1]这个状态,以服从pi[t - 1]为均值,方差1的正态分布的概率转移到了pi_star状态。
pi_star = norm.rvs(loc=pi[t - 1], scale=1, size=1, random_state=None)
# 求alpha。这里假设Q(i,j)与Q(j,i)相同,实际并不相同。不过不影响理解。
# norm_dist_prob(pi_star)即pi_star状态在平稳分布中的概率。只不过我们这里为了说明问题,将平稳分布的
# pdf写出来了,在实际中是不知道的。
alpha = min(1, (norm_dist_prob(pi_star) / norm_dist_prob(pi[t - 1])))
u = random.uniform(0, 1)
if u < alpha:
# 如果均匀采样的值小于alpha,那么就接受这个采样值pi_star。
pi[t] = pi_star
else:
# 如果均匀采样的值大于alpha,即不接受转移,这里的操作是:将t轮的值直接给了t+1轮。
pi[t] = pi[t - 1]
plt.scatter(pi, norm.pdf(pi, loc=3, scale=2))
num_bins = 50
plt.hist(pi, num_bins, normed=1, facecolor='red', alpha=0.7)
plt.show()
注:
① 求解
α
\alpha
α利用的是
α
(
i
,
j
)
=
m
i
n
{
π
(
j
)
π
(
i
)
,
1
}
\alpha(i,j) = min\{ \frac{\pi(j)}{\pi(i)},1\}
α(i,j)=min{π(i)π(j),1},即假设
Q
(
i
,
j
)
=
Q
(
j
,
i
)
Q(i,j) = Q(j,i)
Q(i,j)=Q(j,i)。
②代码中没有体现出 n 1 n_1 n1。因为在实际工程应用中,我们一般不判断阈值,都是默认采样到某一个数量,认为就是我们需要的样本。如果计算量大的时候,这个“数量阈值”作用就不大了。
输出的图中可以看到采样值的分布与真实的分布之间的关系如下,采样集还是比较拟合对应分布的。
**总结一下:**在实际应用MH采样时,我们设定一个马尔科夫链状态转移矩阵 Q Q Q, Q Q Q一般是服从正态分布的,这是因为正态分布的性质较好,它的边缘分布条件分布也是正态分布。另外,我们还需要一个不知道其分布的数据集,我们就是要抽样服从这个数据集分布的样本。
五、M-H采样总结
M-H采样完整解决了使用蒙特卡罗方法需要的任意概率分布样本集的问题,因此在实际生产环境得到了广泛的应用。
但是在大数据时代,M-H采样面临着两大难题:
1) 我们的数据特征非常的多,M-H采样由于接受率计算式 π ( j ) Q ( j , i ) π ( i ) Q ( i , j ) \frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)} π(i)Q(i,j)π(j)Q(j,i)的存在,在高维时需要的计算时间非常的可观,算法效率很低。同时 α ( i , j ) \alpha(i,j) α(i,j)一般小于1,有时候辛苦计算出来却被拒绝了。能不能做到不拒绝转移呢?
2) 由于特征维度大,很多时候我们甚至很难求出目标的各特征维度联合分布,但是可以方便求出各个特征之间的条件概率分布。这时候我们能不能只有各维度之间条件概率分布的情况下方便的采样呢?
Gibbs采样解决了上面两个问题,因此在大数据时代,MCMC采样基本是Gibbs采样的天下。我们在《
MCMC(三):Gibbs采样》来讨论Gibbs采样。