浅谈马尔科夫链蒙特卡罗方法（MCMC）算法的理解

1.解决的问题

计算机怎么在任意给定的概率分布P上采样？首先可以想到把它拆成两步：

（1）首先等概率的从采样区间里取一个待定样本x，并得到它的概率为p(x)

（2）然后在均匀分布U[0,1]上取一个值，如果这个值高于了样本出现的概率p(x)，那么就舍弃掉这个待定样本，如果低于或等于这个值就保留。

不断重复这个步骤，我们就能得到满足概率分布P的样本X。

这个方法看似有效，但存在一个致命的问题，那就是慢，因为概率p的均值，随着采样的区域增大而减小【如果P是多维概率分布，那么随着维度不断增大，那么这个值将趋近于0，也就是维度灾难】，那么就意味着这个样本被保留下来的概率非常小。

那么怎么改进？

我们把上述的第（2）步改一下，不再用均匀U[0,1]进行采样，而是用U[0,max(p)]，而选取的规则变为，如果从均匀分布U[0，max(p)]里取得的值u高于p(x)，那么就保留该样本。

这样能够保证，至少有一个可能出现的样本在概率区间接受率为1。

这种改进比第一种的命中率要高上不少，如果p是均匀分布，那么采样命中率将是100%，当然如果采样是均匀分布，那么我直接取就好了，没必要多这么一步。

那么在此之上还有什么改进呢？

你自然能想到，我们不再用均匀分布U来计算样本是否可以保留，我们找一个好实现的概率分布Q，这个Q要和分布P相似，我们放缩该概率密度函数，使得min(q)>max(p)。

这样一来，命中率就又提高了，而且可以预想到，如果概率分布q=p那么采样命中率就会是100%，但这自然不可能。

那么有没有一种方法可以让采样的命中率是100%呢？

有，就是MCMC，在MCMC中，所有样本都是有效采样。而原理则是用频率取代概率。

2.MCMC

假设我们手里有个小球，每隔一秒就会变一种颜色（包括当前颜色），一段时间后，我们统计每种颜色的时间，就可以得到，各颜色出现的概率。

而这就是MCMC的原理，我们要做的是借助马尔可夫链，实现这样一个小球。

我们给小球输入一个指令，给定它每种颜色出现的概率，也就是P。

这时候小球内部有个状态机会来实现这个概率，首先会按照颜色的数量建立状态，同时每种状态之间会有线路连接，且每种状态之间都是相互可达的，即连通。

这时候有钢珠在每种状态里游走，钢珠到了哪种状态，哪种状态代表的灯就会亮起。

好了，如何让小球的灯能够以某种稳定的频率让灯亮起。

研究发现，只要让两两状态之间转换的概率与对方状态的概率乘积相等。

即P(A)Q(A->B)=P(B)Q(B->A)，这就是马尔科夫链细致平衡条件，A和B代表任意两个状态。Q(A-B)是状态A到状态B的转移概率。

这时候，我们只要跟随钢珠在各状态之间的游走，记录下状态值即样本值，就能完成采样。

具体步骤如下，

（1）记录钢珠初始值状态X1。

（2）随机选择一个目标状态。

（3）查看抵达目标状态的转移概率。

（4）从U[0,1]中采一个值，如果值大于转移概率，那么状态不变，如果小于或等于转移概率，那么状态转化为目标状态。

重复如上步骤，就能得到一个服从P分布的样本了。

但这个算法虽然没有无效采样，但状态A转移到状态的转移概率过低，那么就会长时间卡在这一种状态。那么有没有方法改进呢。

当然有，那就是等比例缩放平衡方程两边的转移概率，P(A)Q(A->B)=P(B)Q(B->A)。

这里有个限制，那就是不能破坏转移概率小于或等于1概率特性，同时也不能破坏平衡方程本身。

那么我们能找到的一个缩放系数就是MAX[P(B->A),P(A->B)]，两边同时除以该系数。

则有P(A)MIN[1,Q(A->B)/Q(B->A)] = P(B)MIN[1,Q(A->B)/Q(B->A)]

这样就能保证至少有一边可以必定转化为另一边。

到这里MCMC就算比较完整了。

问题就是Q怎么选。

这里有两种方案，一个就是Q任取，只要满足转移矩阵的要求，然后补一个系数配平方程

即: P(A)Q(A->B)*[P(B)Q(B->A)]=P(B)Q(B->A)*[P(A)Q(A->B)]，

缩放后有：

P(A)Q(B->A)*MIN[1,P(B)Q(B->A)/P(B)Q(B->A)] = P(B)Q(B->A)*MIN[1,P(A)Q(A->B)/P(B)Q(B->A)]

另一种那就是选Q(A->B) = P(B)，Q(B->A)=P(A)。

前者是HM的思想，配平的系数的使用方法类比前一节的转移接受率。

后者是Gibbs的思想。

具体实现有很多细节，可以看下面推荐的教学视频，这里推荐先看P6

 机器学习-白板推导系列-蒙特卡洛方法5（Monte Carlo Method）- 再回首_哔哩哔哩_bilibili