在介绍《Fast and Provably Good Seedings for K-means》时，作者使用MCMC采样来近似D2−sampling 的过程，下面我们来介绍一下MCMC算法。



Monte Carlo Approach

在介绍MCMC算法之前，我们先来看蒙特卡洛随机模拟算法。

假设我们需要求解下面这个积分问题： 
m=∫baf(x)d(x)
而这个f(x) 的积分形式比较难求，我们可以通过数值模拟解法近似求解。常用的方法是蒙特卡洛算法： 
m=∫baf(x)q(x)q(x)d(x)
式中的q(x) 可以看作 x 在区间上的分布，当我们使用 q(x) 采集足够的样本后，便可以使用均值来近似积分： 
m=1n∑i=1Nf(xi)q(xi)
以上便是蒙特卡洛算法的核心思想，具体操作时就是需要考虑如何从给定的概率分布 (x) 中采集足够的样本。下面我们来介绍两种简单的采样算法。



采样算法

CDF

假设我们想以 p(x) 的概率采集一些数据，


  
我们可以先求出它的累积概率函数: 
F(x)=∫bap(x)d(x)   s.t.a<x<b


  2.然后使用均匀分布产生一个样本点 u
u∼u(a,b)


  3.最后我们求解出反函数F−1(u)
x=F−1(u)


  4.x 便可以看作是按照 p(x) 采样得到的样本点
  

我们可以看一个正态分布的例子：


  

  
  
左图是概率密度函数PDF-p(x)，右图是对应的累计概率密度函数CDF-F(x). 

  我们先随机生成一个样本点 u=0.8413， F−1(u)=1，所以1就是相应的样本点。 

  直觉上看，我们随机选择(0,1)之间的点，大部分点都i会落在区间(0.2,0.8) 之间，而与之对应的 F−1(x) 则大部分落在(−1,1)之间，这也对应了正态分布在(−1,1)之间概率密度较大的事实。


这种方法只适用于简单的分布。

Markov Chain

对于某一些马尔科夫链，它能够收敛到一个平稳分布π(x)，即： 
π=πP


当马尔可夫链在第 n 个状态收敛之后，从 n+1 往后的状态都可以看作是平稳分布 π 的采样点。 

所以为了构造一个以 π 为平稳分布的马尔科夫链，我们需要满足细致平衡条件:



π(i)P(i,j)=π(j)P(j,i)
Metropolis Sampler



M算法要求提案概率 q(x) 是对称的。 
Metropolis-Hastings Sampler



MH算法加入了 α ，不要求提案概率是对称的，也能满足细致平衡条件。

这两个算法的区别在于提案概率是否是对称的。

统计学习-MCMC采样与MH采样简介
初探MCMC采样方法和MH采样方法

    在之前介绍的文章中，我们介绍了基本的采样方法。在了解简单的采样方法后，本次文章将会介绍MCMC采样方法与优化的MH采样方法，这两种方法对于一些高维度特征的数据可以进行比较合理的采样。

    MCMC采样方法其实就是蒙特卡洛马尔科夫链的采样方法，因为这两个术语的英文简写都是MC，所以我们使用MCMC采样来代替。MH采样方法的英文是Metropolis-Hasting，它是一种基于MCMC采样的优化方法，在实际机器学习问题中使用较多。下面我们开始介绍这两种方法：

    首先是MCMC方法，该方法里的蒙特卡洛方法是一种计算数值积分的方法，将一些很难从解析形式求解的定积分转为基于概率分布的级数求和的方法进行逼近求解。但是要使用蒙特卡洛方法最重要的是需要事先明确概率分布的函数，如果是简单的分布形式比如F分布、Gamma分布等，可以通过与均匀分布进行映射变换得到；而对于复杂分布或者不常见分布，则需要进行“接受-拒绝采样”来进行逼近原始的复杂分布，这种思路就是用一个简单的分布形式比如正态分布去完全覆盖原有的复杂分布，也就是下图的两个分布形式，如果每次抽样落在灰色区域的时候就拒绝抽样，反之则接受抽样。



    其实要介绍的是MCMC方法里的马尔科夫链方法，这种马尔可夫链假设的是当前状态只和之前的一个状态有关系存在，和之前的历史状态都无关。那么这种非周期的马尔可夫链就有一个比较好的性质，也就是会收敛到一个稳态的分布。使用马尔可夫链方法进行抽样的时候，只需要明确状态转移矩阵和初始的状态，就可以得出后续的样本点，从而得到所需的采样样本。对于马尔可夫链其实还有一种非常重要的细致平稳条件，在所有的状态里任意取i、j两个状态，都满足G(i)P(i,j)=G(j)P(j,i)  ，而MCMC采样和MH采样的基本原理也是根据这个性质延伸得到的。

    对于MCMC采样方法，我们可以在细致平稳条件加上概率b(i,j)得到，也就是将原来的等式约束进行放松，得到G(i)P(i,j)b(i,j)=G(j)P(j,i)b(j,i)  ，所以整个MCMC采样的步骤如下所示：

    1）初始化相关的采样参数，比如需要的样本数N，状态转移矩阵T，平稳分布G，转移状态次数的阈值U等参数

    2）从一个任意的简单的分布进行采样，得到初始状态值x0

    3）对于当前状态t=0到N+U-1的状态之中，重复下述的4-5-6步骤：

    4）从条件概率分布T(x|x(t))中进行采样得到样本x(k)

    5）从均匀分布进行采样，得到样本v

    6）使用之前的宽松的细致平稳条件进行判断，即样本v是否满足小于G(x(k))T(x(k),x(t)) ，如果满足的话就接受转移，并且让x(t+1)=x(k) ，否则就不接受转移，让x(t+1)=x(t)

    通过上述步骤后最终得到的就是基于平稳分布的样本集。但是由于接受的转移概率比较小，所以这种MCMC方法很难得到真正的样本，因为很多样本都会被拒绝转移，整体的采样效率较低。

    对于MH方法，其实步骤和上述的MCMC采样方法相似，它只是将阈值放大了5倍，对于阈值b(i,j)的公式进行修正：b(i,j)=min(G(j)T(j,i)/G(i)T(i,j),1) 从而得到MH采样方法，其他步骤和MCMC方法完全一样。但是需要注意的是，MH采样和MCMC方法一样，对于特征较多的情况也是比较难使用的。

    总的来说，MCMC采样方法和MH采样方法都是比较合理的采样方法，但是MH采样是MCMC采样方法的简化版，所以在实际问题中使用MH采样的比较多。MCMC采样在数学上是比较严谨和稳定的，但是满足它采样接受率的样本比较少，也就是说总是会采样失败，这种时候对原来的方法进行改进以后也就是MH方法，可以有效避免上述问题。MCMC方法和MH采样是在数据特征不多的时候，被广泛使用的采样方法，当维度变得很多的时候会使用其他的采样方法，这个在后续的文章会进行介绍。而对于初学者则需要了解这类采样方法的原理，方便在实际问题中更好地进行使用。

　　　MCMC(一)蒙特卡罗方法

　　　　MCMC(二)马尔科夫链

　　　　MCMC(三)MCMC采样和M-H采样

　　　　MCMC(四)Gibbs采样

 

　　　　在MCMC(二)马尔科夫链中我们讲到给定一个概率平稳分布ππ, 很难直接找到对应的马尔科夫链状态转移矩阵PP。而只要解决这个问题，我们就可以找到一种通用的概率分布采样方法，进而用于蒙特卡罗模拟。本篇我们就讨论解决这个问题的办法：MCMC采样和它的易用版M-H采样。

1. 马尔科夫链的细致平稳条件

　　　　在解决从平稳分布ππ, 找到对应的马尔科夫链状态转移矩阵PP之前，我们还需要先看看马尔科夫链的细致平稳条件。定义如下：

　　　　如果非周期马尔科夫链的状态转移矩阵PP和概率分布π(x)π(x)对于所有的i,ji,j满足：

π(i)P(i,j)=π(j)P(j,i)π(i)P(i,j)=π(j)P(j,i)

　　　　则称概率分布π(x)π(x)是状态转移矩阵PP的平稳分布。

　　　　证明很简单,由细致平稳条件有：

∑i=1∞π(i)P(i,j)=∑i=1∞π(j)P(j,i)=π(j)∑i=1∞P(j,i)=π(j)∑i=1∞π(i)P(i,j)=∑i=1∞π(j)P(j,i)=π(j)∑i=1∞P(j,i)=π(j)

　　　　将上式用矩阵表示即为：

πP=ππP=π

　　　　即满足马尔可夫链的收敛性质。也就是说，只要我们找到了可以使概率分布π(x)π(x)满足细致平稳分布的矩阵PP即可。这给了我们寻找从平稳分布ππ, 找到对应的马尔科夫链状态转移矩阵PP的新思路。

　　　　不过不幸的是，仅仅从细致平稳条件还是很难找到合适的矩阵PP。比如我们的目标平稳分布是π(x)π(x),随机找一个马尔科夫链状态转移矩阵QQ,它是很难满足细致平稳条件的，即：

π(i)Q(i,j)≠π(j)Q(j,i)π(i)Q(i,j)≠π(j)Q(j,i)

　　　　那么如何使这个等式满足呢？下面我们来看MCMC采样如何解决这个问题。

2. MCMC采样

　　　　由于一般情况下，目标平稳分布π(x)π(x)和某一个马尔科夫链状态转移矩阵QQ不满足细致平稳条件，即

π(i)Q(i,j)≠π(j)Q(j,i)π(i)Q(i,j)≠π(j)Q(j,i)

　　　　我们可以对上式做一个改造，使细致平稳条件成立。方法是引入一个α(i,j)α(i,j),使上式可以取等号，即：

π(i)Q(i,j)α(i,j)=π(j)Q(j,i)α(j,i)π(i)Q(i,j)α(i,j)=π(j)Q(j,i)α(j,i)

　　　　问题是什么样的α(i,j)α(i,j)可以使等式成立呢？其实很简单，只要满足下两式即可：

α(i,j)=π(j)Q(j,i)α(i,j)=π(j)Q(j,i)

α(j,i)=π(i)Q(i,j)α(j,i)=π(i)Q(i,j)

　　　　这样，我们就得到了我们的分布π(x)π(x)对应的马尔科夫链状态转移矩阵PP，满足：

P(i,j)=Q(i,j)α(i,j)P(i,j)=Q(i,j)α(i,j)

　　　　也就是说，我们的目标矩阵PP可以通过任意一个马尔科夫链状态转移矩阵QQ乘以α(i,j)α(i,j)得到。α(i,j)α(i,j)我们有一般称之为接受率。取值在[0,1]之间，可以理解为一个概率值。即目标矩阵PP可以通过任意一个马尔科夫链状态转移矩阵QQ以一定的接受率获得。这个很像我们在MCMC(一)蒙特卡罗方法第4节讲到的接受-拒绝采样，那里是以一个常用分布通过一定的接受-拒绝概率得到一个非常见分布，这里是以一个常见的马尔科夫链状态转移矩阵QQ通过一定的接受-拒绝概率得到目标转移矩阵PP,两者的解决问题思路是类似的。

　　　　好了，现在我们来总结下MCMC的采样过程。

　　　　1）输入我们任意选定的马尔科夫链状态转移矩阵QQ，平稳分布π(x)π(x)，设定状态转移次数阈值n1n1，需要的样本个数n2n2

　　　　2）从任意简单概率分布采样得到初始状态值x0x0

　　　　3）for t=0t=0 to n1+n2−1n1+n2−1: 

　　　　　　a) 从条件概率分布Q(x|xt)Q(x|xt)中采样得到样本x∗x∗

　　　　　　b) 从均匀分布采样u∼uniform[0,1]u∼uniform[0,1]

　　　　　　c) 如果u<α(xt,x∗)=π(x∗)Q(x∗,xt)u<α(xt,x∗)=π(x∗)Q(x∗,xt), 则接受转移xt→x∗xt→x∗，即xt+1=x∗xt+1=x∗

　　　　　　d) 否则不接受转移，即xt+1=xtxt+1=xt

　　　　样本集(xn1,xn1+1,...,xn1+n2−1)(xn1,xn1+1,...,xn1+n2−1)即为我们需要的平稳分布对应的样本集。

　　　　上面这个过程基本上就是MCMC采样的完整采样理论了，但是这个采样算法还是比较难在实际中应用，为什么呢？问题在上面第三步的c步骤，接受率这儿。由于α(xt,x∗)α(xt,x∗)可能非常的小，比如0.1，导致我们大部分的采样值都被拒绝转移，采样效率很低。有可能我们采样了上百万次马尔可夫链还没有收敛，也就是上面这个n1n1要非常非常的大，这让人难以接受，怎么办呢？这时就轮到我们的M-H采样出场了。

3. 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)π(i)Q(i,j)α(i,j)=π(j)Q(j,i)α(j,i)

　　　　我们采样效率低的原因是α(i,j)α(i,j)太小了，比如为0.1，而α(j,i)α(j,i)为0.2。即：

π(i)Q(i,j)×0.1=π(j)Q(j,i)×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π(i)Q(i,j)×0.5=π(j)Q(j,i)×1

　　　　这样我们的接受率可以做如下改进，即：

α(i,j)=min{π(j)Q(j,i)π(i)Q(i,j),1}α(i,j)=min{π(j)Q(j,i)π(i)Q(i,j),1}

　　　　通过这个微小的改造，我们就得到了可以在实际应用中使用的M-H采样算法过程如下：

　　　　1）输入我们任意选定的马尔科夫链状态转移矩阵QQ，平稳分布π(x)π(x)，设定状态转移次数阈值n1n1，需要的样本个数n2n2

　　　　2）从任意简单概率分布采样得到初始状态值x0x0

　　　　3）for t=0t=0 to n1+n2−1n1+n2−1: 

　　　　　　a) 从条件概率分布Q(x|xt)Q(x|xt)中采样得到样本x∗x∗

　　　　　　b) 从均匀分布采样u∼uniform[0,1]u∼uniform[0,1]

　　　　　　c) 如果u<α(xt,x∗)=min{π(j)Q(j,i)π(i)Q(i,j),1}u<α(xt,x∗)=min{π(j)Q(j,i)π(i)Q(i,j),1}, 则接受转移xt→x∗xt→x∗，即xt+1=x∗xt+1=x∗

　　　　　　d) 否则不接受转移，即xt+1=xtxt+1=xt

　　　　样本集(xn1,xn1+1,...,xn1+n2−1)(xn1,xn1+1,...,xn1+n2−1)即为我们需要的平稳分布对应的样本集。

　　　　很多时候，我们选择的马尔科夫链状态转移矩阵QQ如果是对称的，即满足Q(i,j)=Q(j,i)Q(i,j)=Q(j,i),这时我们的接受率可以进一步简化为：

α(i,j)=min{π(j)π(i),1}α(i,j)=min{π(j)π(i),1}

4. M-H采样实例

　　　　为了更容易理解，这里给出一个M-H采样的实例。

　　　　完整代码参见我的github: https://github.com/ljpzzz/machinelearning/blob/master/mathematics/mcmc_3_4.ipynb

　　　　在例子里，我们的目标平稳分布是一个均值3，标准差2的正态分布，而选择的马尔可夫链状态转移矩阵Q(i,j)Q(i,j)的条件转移概率是以ii为均值,方差1的正态分布在位置jj的值。这个例子仅仅用来让大家加深对M-H采样过程的理解。毕竟一个普通的一维正态分布用不着去用M-H采样来获得样本。

　　　　代码如下：



import random
import math
from scipy.stats import norm
import matplotlib.pyplot as plt
%matplotlib inline

def norm_dist_prob(theta):
    y = norm.pdf(theta, loc=3, scale=2)
    return y

T = 5000
pi = [0 for i in range(T)]
sigma = 1
t = 0
while t < T-1:
    t = t + 1
    pi_star = norm.rvs(loc=pi[t - 1], scale=sigma, size=1, random_state=None)
    alpha = min(1, (norm_dist_prob(pi_star[0]) / norm_dist_prob(pi[t - 1])))

    u = random.uniform(0, 1)
    if u < alpha:
        pi[t] = pi_star[0]
    else:
        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()



　　　　输出的图中可以看到采样值的分布与真实的分布之间的关系如下，采样集还是比较拟合对应分布的。



 

5. M-H采样总结

　　　　M-H采样完整解决了使用蒙特卡罗方法需要的任意概率分布样本集的问题，因此在实际生产环境得到了广泛的应用。

　　　　但是在大数据时代，M-H采样面临着两大难题：

　　　　1） 我们的数据特征非常的多，M-H采样由于接受率计算式π(j)Q(j,i)π(i)Q(i,j)π(j)Q(j,i)π(i)Q(i,j)的存在，在高维时需要的计算时间非常的可观，算法效率很低。同时α(i,j)α(i,j)一般小于1，有时候辛苦计算出来却被拒绝了。能不能做到不拒绝转移呢？

　　　　2） 由于特征维度大，很多时候我们甚至很难求出目标的各特征维度联合分布，但是可以方便求出各个特征之间的条件概率分布。这时候我们能不能只有各维度之间条件概率分布的情况下方便的采样呢？

　　　　Gibbs采样解决了上面两个问题，因此在大数据时代，MCMC采样基本是Gibbs采样的天下，下一篇我们就来讨论Gibbs采样。

MCMC采样算法

完整的MCMC采样算法已经有很多博主发布了，这里就不再重复了。主要想分享一下在看其他博主写的MCMC采样算法时，不太理解的地方。

MCMC采样关键问题在于如何构建转移矩阵，使得平稳分布恰好是p(x)。主要使用细致平稳条件。



细致平稳条件

如果非周期马氏链的转移矩阵P和分布π(x)满足： 

π(i)Pij=π(j)Pji  for all i,j 

则π(x)是马尔可夫链的平稳分布，上式称为细致平稳条件。

马氏链的一个例子： 

社会学家经常把人按其经济状况分成3类：下层(lower-class)、中层(middle-class)、上层(upper-class)，我们用1,2,3 分别代表这三个阶层。社会学家们发现决定一个人的收入阶层的最重要的因素就是其父母的收入阶层。如果一个人的收入属于下层类别，那么他的孩子属于下层收入的概率是 0.65, 属于中层收入的概率是 0.28, 属于上层收入的概率是 0.07。事实上，从父代到子代，收入阶层的变化的转移概率如下 




假设初始概率分布为π0=[0.21,0.68,0.11]，则我们可以计算前n代人的分布状况如下： 


从第7代人开始，这个分布开始稳定不变。

回到细致平稳分布，当达到马氏链的平稳分布的时候，π=(0.286，0.489，0.225)，由细致平稳可知满足π(i)Pij = π(j)Pji。 

如： π(1)P12=π(2)P21—–>0.286*0.28 =0.08008 

                                              0.489*0.15 =0.07335  近似相等 

所以当π(i)Pij=π(j)Pji时，π(x)是马氏链的平稳分布。



MCMC算法

假设我们已经有一个转移矩阵Q(对应元素为q(i,j)), 得到了如下的用于采样概率分布p(x)的算法。 


读Rickjin的LDA数学八卦的时候特别不理解，为什么α就是接受率。 

先谈论一个简单的采样，假设采样掷一枚硬币是正面还是反面，从均匀分布中采样u，当u>0.5时候接受采样（正面），否则拒绝采样。

下面说说我对MCMC算法的理解： 

假设对一篇文章的4个主题分布进行采样：设置初始状态X0=Topic1。 

假设时刻t=1,从转移矩阵Q中的q(x|T1)分布采样了T2。 

对于随机游走，此时以概率p(T1)q(T2|T1)从Topic1转移到Topic2，采样T2，没有接受率之说。 

而此时我们需要采样的是马氏链细致平稳分布π(x)/p(x)的样本，必须满足π(i)Pij=π(j)Pji。为了保证p(T1)q(T2|T1)=p(T2)q(T1|T2)，使得样本来自p(x)概率分布,只有当u小于p(T2)q(T1|T2)时，接受X1=Topic2。

MCMC理解（网上看到的另一种解释）




Gibbs sampling

•   什么是Gibbs Sampling 

Gibbs Sampling是MCMC算法中的一种，用来构造多变量概率分布的随机样本，比如构造两个或多个变量的联合分布，求积分，期望。 

•   为什么需要Gibbs Sampling 

积分、期望、联合分布计算，通常情况下当前面三个问题是NP问题时才需要Gibbs Sampling，否则直接计算就可以了。补充一句Gibbs Sampling只是（也只能）到近似解。 

•   应用场景 

a、积分，期望，样本概率很难计算出来；b、条件概率很容易计算。具体例子：受限玻尔兹曼机（RBM）的训练，贝叶斯网络，LDA都用到Gibbs Sampling。 

•   为什么Gibbs Sampling有效 

当Gibbs Sapling算法执行多次之后，产生的样本服从真实样本的分布，即相当于直接从联合分布中采样。