在前面《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 pltdef 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_starelse:# 如果均匀采样的值大于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采样。

参考文献

【1】MCMC(三)MCMC采样和M-H采样

【2】随机采样方法整理与讲解（MCMC、Gibbs Sampling等）

MCMC(二)：MCMC采样和M-H采样相关推荐

漫谈MCMC与Gibbs采样（一）—— 采样背后的逻辑
采样背后的逻辑我们为什么需要采样我们先来思考一个最基本的问题,我们为什么需要采样?之所以讨论这一点,是因为深刻理解我们要解决的问题,可以帮助我们更好地掌握解决问题的方法.尤其是当我们在复杂的模型中 ...
数字图像处理（二）灰度、量化、采样
转自https://blog.csdn.net/eastmount/article/details/46010637 本文主要讲述基于VC++6.0 MFC图像处理的应用知识,主要结合自己大三所学课程 ...
深度学习 --- 受限玻尔兹曼机RBM（直接采样、接受-拒绝采样、重要性采样详解）
在讲解MCMC和Gibbs采样之前,大家需要理解统计学中的采样,什么是采样?为什么要采样?采样有什么用?大家需要深入理解采样的原理,深入理解的好处不仅容易理解下面的MCMC和Gibbs采样,也更容易掌 ...
重要性采样和多重重要性采样在路径追踪中的应用
重要性采样和多重重要性采样在路径追踪中的应用 1 蒙特卡洛路径追踪简要回顾 1.1 算法主要流程 1.2 半球面均匀采样方法 2 重要性采样的运用 2.1 简单例子与基本概念 2.2 路径追踪中的重要 ...
粒子滤波 particle filter —从贝叶斯滤波到粒子滤波—Part-III（重要性采样序贯重要性采样SIS）
粒子滤波 particle filter -从贝叶斯滤波到粒子滤波-Part-III(重要性采样&序贯重要性采样SIS) 原创不易,路过的各位大佬请点个赞机动目标跟踪/非线性滤波/传感器融合 ...
为什么先编码再解码？即先降采样，然后上采样
目录下采样原理: 回答一: 回答二: 回答三: 下采样原理: 对于一幅图像I尺寸为M*N,对其进行s倍下采样,即得到(M/s)*(N/s)尺寸的得分辨率图像,当然s应该是M和N的公约数才行,如果考虑 ...
[Matlab-4]信号的采样与恢复（采样定理）
[Matlab-4]信号的采样与恢复(采样定理) 采样定理(Nyquist Sampling Theorem) 移动平均滤波器(Moving-average filter) 巴特沃斯滤波器(Butte ...
pcm 降采样_图像降采样和升采样
转自:http://www.lofter.com/postentry?from=search&permalink=1cb3111d_6ee9587 1.先说说这两个词的概念: 降采样,即是采样 ...
DataScience：对严重不均衡数据集进行多种采样策略(随机过抽样、SMOTE过采样、SMOTETomek综合采样、改变样本权重等)简介、经验总结之详细攻略
DataScience:对严重不均衡数据集进行多种采样策略(随机过抽样.SMOTE过采样.SMOTETomek综合采样.改变样本权重等)简介.经验总结之详细攻略目录

MCMC(二)：MCMC采样和M-H采样