Metropolis-Hasting抽样算法

1. 随机模拟的基本思想

假设我们有一个矩形区域RRR，面积为S0S_0S0。在此区域中，有一个不规则区域MMM，其面积SSS待求。

方法1：把不规则区域MMM划分为多个小的规则区域，由这些规则区域的面积总和S′S'S′近似。

方法2：抓一把黄豆，均匀铺在RRR中，再统计落在不规则区域MMM中的黄豆比例。该比例可近似SS0\frac {S} {S_0}S0S。
方法2采用的是抽样（采样）解决问题的思想，即随机模拟。

因此，随机模拟的关键，在于如何将实际复杂问题，转化为抽样可以解决的问题。当然，这非常灵活，也考验我们的创造力。

我们接下来的核心问题，是利用计算机可直接实现的均匀抽样，借助随机模拟的方法，实现满足特定概率分布p(x)p(x)p(x)的抽样。

2. 拒绝抽样

当我们的概率分布很简单时，如[0,1][0,1][0,1]上的均匀分布，抽样是很好实现的。现成的随机算法非常多。

但是，对于其他更复杂的分布，如果我们还想借助简单的均匀随机抽样，就需要换换思路了。

例如上图中的p(z)p(z)p(z)，是较为复杂的概率密度函数。

我们先构造比较简单的、容易抽样的高斯分布q(z)q(z)q(z)，乘一个常数kkk，使之满足：
kq(z)⩾p(z),∀zkq(z) \geqslant p(z), \forall z kq(z)⩾p(z),∀z

拒绝抽样的方法是：

对高斯分布的样本进行采样，得到一个z0z_0z0，如图；
在[0,kq(z0)][0,kq(z_0)][0,kq(z0)]上均匀采样，得到u0u_0u0；
如果u0>p(z0)u_0>p(z_0)u0>p(z0)，就拒绝该采样结果；反之接受；
重复直到接受足够多的样本。

我们理解一下为什么拒绝抽样是可行的，这对我们之后运用随机模拟思想非常关键！

该简单分布与高斯分布的趋势大致相同，中间大两头小，因此取到中间的zzz的概率比较大；
通过u0u_0u0进一步筛选zzz，p(z)p(z)p(z)较低时拒绝率较大，反之接受率较高，大致符合p(z)p(z)p(z)分布要求。

一句话，通过简单抽样和设置接受率，“强迫”结果趋于复杂分布。

这样，我们就通过简单的高斯分布采样，以及计算机可直接实现的均匀采样，间接实现了复杂的p(z)p(z)p(z)采样。

然而，在高维情况下，该方法的劣势很明显：

合适的简单分布q(z)q(z)q(z)很难找到。
合适的kkk值很难找到。

这两个问题会导致拒绝率极高，计算效率很低。其根本缺陷在于：样本之间是独立无关的。

3. Metropolis-Hastings抽样

3.1 引入思想

结合马氏链知识（参考信息论相关书籍），我们知道：
对于遍历的马尔可夫链，当其转移次数足够多时，将会进入稳态分布π(x)\pi (x)π(x)，即各状态的出现概率趋于稳定。

进一步，假设变量XXX所有可能的取值，构成某遍历马氏链的状态空间。
则无论从什么状态出发，只要转移步数足够大，该马氏链将趋于稳定，即逐渐开始依次输出符合p(x)p(x)p(x)分布的状态XiX_iXi。
这样，通过收集马氏链收敛后产生的XiX_iXi，我们就得到了符合分布p(x)p(x)p(x)的样本。
此时稳态分布即π(x)=p(x)\pi(x)=p(x)π(x)=p(x)。

比如，我们希望抽样样本满足指数分布。此时，整数0到正无穷都是状态，但整数0的出现概率最大，随着数增大，出现概率越来越小。我们构造一个稳定输出服从指数分布状态的马氏链，则得到的样本服从指数分布。

这个绝妙的想法在1953年由Metropolis提出。为了研究粒子系统的平稳性质，Metropolis 考虑了物理学中常见的波尔兹曼分布的采样问题，首次提出了基于马氏链的蒙特卡罗方法（Metropolis算法），并在最早的计算机上编程实现。

Metropolis 算法是首个普适的采样方法，并启发了一系列 MCMC方法，所以人们把它视为随机模拟技术腾飞的起点。 Metropolis的这篇论文被收录在《统计学中的重大突破》中， Metropolis算法也被选为二十世纪的十个最重要的算法之一。

我们接下来介绍的MH算法是Metropolis 算法的一个常用的改进变种。
由于马氏链的收敛性质主要由转移矩阵PPP决定, 所以基于马氏链做采样的关键问题是：如何构造转移矩阵PPP，使平稳分布恰好是我们要的分布p(x)p(x)p(x)。

3.2 理论基础：细致平稳条件

定理——细致平稳条件：

若非周期马氏链的转移矩阵PPP和分布π(x)\pi(x)π(x)满足：
π(i)Pij=π(j)Pji,∀i,j\pi(i)P_{ij}=\pi(j)P_{ji}, \forall i,j π(i)Pij=π(j)Pji,∀i,j
则分布π(x)\pi(x)π(x)是马氏链的平稳分布。

证明：
∑i=1∞π(i)Pij=∑i=1∞π(j)Pji=π(j)∑i=1∞Pji=π(j),∀j\sum_{i=1}^\infty \pi(i)P_{ij} = \sum_{i=1}^\infty \pi(j)P_{ji} = \pi(j)\sum_{i=1}^\infty P_{ji} = \pi(j), \forall j i=1∑∞π(i)Pij=i=1∑∞π(j)Pji=π(j)i=1∑∞Pji=π(j),∀j
这说明：
π(x)P=π(x)\pi(x)P = \pi(x) π(x)P=π(x)
因此π(x)\pi(x)π(x)是平稳分布。

物理意义：

对于任意两个状态i,ji,ji,j，从状态iii流失到状态jjj的概率质量，等于反向流动的概率质量，因此是状态概率是稳定的。

3.3 M-H算法实现

假设我们已经有了一个转移矩阵QQQ。一般情况下，该转移矩阵不满足细致平稳条件：
P(i)Qij≠P(j)QjiP(i)Q_{ij} \not= P(j)Q_{ji} P(i)Qij=P(j)Qji
我们引入接受率α(i,j)\alpha (i,j)α(i,j)，满足：
P(i)Qijα(i,j)=P(j)Qjiα(j,i)P(i)Q_{ij} \alpha (i,j) = P(j)Q_{ji} \alpha (j,i) P(i)Qijα(i,j)=P(j)Qjiα(j,i)
其中：

P(i)P(i)P(i)是状态iii出现的概率，是假设马氏链满足复杂分布时输出状态的概率。
QQQ是初始转移概率矩阵，满足任意一个给定的简单分布，便于抽样即可。

显然，接受率最简单的构造方法是：
α(i,j)=P(j)Qji\alpha(i,j)=P(j)Q_{ji} α(i,j)=P(j)Qji
α(j,i)=P(i)Qij\alpha(j,i)=P(i)Q_{ij} α(j,i)=P(i)Qij

此时，细致平稳条件成立，该马氏链输出的状态满足稳态分布p(x)p(x)p(x)！

综上，有几个要点必须要理解，有助于我们编程实现：

我们引入接受率，使得该马氏链以转移概率Qijα(i,j)Q_{ij}\alpha(i,j)Qijα(i,j)从状态iii转移到状态jjj。
该转移概率的实现方式是：我们先按QijQ_{ij}Qij生成新状态，再按α(i,j)\alpha(i,j)α(i,j)的概率接受转移结果，乘积即为Qijα(i,j)Q_{ij}\alpha(i,j)Qijα(i,j)转移概率。
一定要理解这个接受率！在拒绝抽样中，如果u0>p(z)u_0>p(z)u0>p(z)，我们会放弃抽样结果，不纳入最终的统计。正是因为如此，抽样才会逼近复杂分布。而这里的接受是“接受转移”，如果u>αu>\alphau>α，意味着t+1t+1t+1时刻状态和ttt时刻相同，并且应该纳入统计而不是放弃结果！！
由于状态出现概率和转移概率满足细致平稳条件，因此状态出现概率即稳态分布概率。长期转移后，我们会得到想要的抽样分布效果。

图解转移过程：

要注意的是，当马氏链处在状态iii时，我们并不知道如何选择下一个状态jjj。
我们在这里采用抽样的方法，并借助接受率，“强迫”转移过程逼近理想的遍历马氏链转移过程。

算法描述如下：

初始化马氏链状态X0=x0X_0=x_0X0=x0
从t=0t=0t=0开始，重复以下步骤：
- 假设该ttt时刻状态为Xt=xtX_t=x_tXt=xt；
- 对Qxt,xQ_{x_t,x}Qxt,x抽样，随机得到一个状态xnextx_{next}xnext；
- 在[0,1][0,1][0,1]上均匀抽样，得到一个uuu；
- 若u<α(xt,xnext)=P(xnext)Qxnext,xtu<\alpha(x_t,x_{next})=P(x_{next})Q_{x_{next},x_t}u<α(xt,xnext)=P(xnext)Qxnext,xt，则接受转移Xt+1=xnextX_{t+1}=x_{next}Xt+1=xnext；
- 否则不转移，即Xt+1=xtX_{t+1}=x_tXt+1=xt。

3.4 算法升级

如果α(xt,xnext)\alpha(x_t,x_{next})α(xt,xnext)太小，会导致转移很难发生，马氏链收敛过慢。

我们回到细致平稳条件：
P(i)Qijα(i,j)=P(j)Qjiα(j,i)P(i)Q_{ij}\alpha(i,j)=P(j)Q_{ji}\alpha(j,i) P(i)Qijα(i,j)=P(j)Qjiα(j,i)
我们希望在保持条件成立的前提下，使接受率尽量增大。
由于接受率的本质是概率，因此α(i,j)\alpha(i,j)α(i,j)和α(j,i)\alpha(j,i)α(j,i)中的较大者不应大于1。那么我们作下述改进即可：
β=P(j)QjiP(i)Qij\beta = \frac {P(j)Q_{ji}} {P(i)Q_{ij}} β=P(i)QijP(j)Qji
α(i,j)=min(β,1)\alpha(i,j)=min(\beta,1) α(i,j)=min(β,1)

解释：

如果α(i,j)\alpha(i,j)α(i,j)更大，那么应设为1。而第一项大于1，因此α(i,j)=1\alpha(i,j)=1α(i,j)=1，结束。
如果α(j,i)\alpha(j,i)α(j,i)更大，那么为了等式成立，α(i,j)\alpha(i,j)α(i,j)必须等于β\betaβ。而β<1\beta<1β<1，得逞～

综上，我们得到了升级版算法：

3.5 仿真实验

仿真目标：

用标准正态分布模拟指数分布（λ=0.1,μ=1λ=10\lambda=0.1,\mu=\frac 1 \lambda=10λ=0.1,μ=λ1=10）。

简化：

由于高斯分布是一个径向基函数（取值仅仅依赖于离原点距离的实值函数），因此该矩阵是转置对称的，Qij=QjiQ_{ij}=Q_{ji}Qij=Qji

代码实现：

#include "cvplot.h"
#include <random>
#include <vector>
#include <algorithm>
using namespace std;double exppdf(double lambda, double x)
{return lambda / exp(lambda * x);
}double exppdf_div(double lambda, double x1, double x2)
{return exp(lambda * (x2 - x1));
}vector<double> bins(vector<double> values, const int nb)
{double min = values[0];double max = values[0];for(double v : values){if (v < min){min = v;}if (v > max){max = v;}}double interval = (max - min) / nb;vector<double> freq(nb + 1, 0);if (interval > 0){for (auto i = 0; i < values.size(); ++i){++freq[(int)(values[i] / interval)];}}return freq;
}int main(void)
{const double lambda = 20;default_random_engine engine;exponential_distribution<double> exp_dist(lambda);uniform_real_distribution<double> unif;size_t i = 0;const int N = 20000;vector<double> ref(N);vector<double> sample(N);while (++i < N){ref[i] = exp_dist(engine);}double current = 1;double u;double alpha;double beta;double next;i = 0;while (i < N){normal_distribution<double> norm_dist(current, 1);next = fabs(norm_dist(engine));beta = exppdf_div(lambda, next, current);alpha = beta < 1 ? beta : 1;u = unif(engine);if (u < alpha){sample[i] = next;}else{sample[i] = current;}current = sample[i];++i;}auto bin1 = bins(ref, 20);auto bin2 = bins(sample,20);cvplot::Series s1("ref", cvplot::chart::Bar);s1.SetRenderColor(cvplot::color::Orange);cvplot::Series s2("sample", cvplot::chart::Bar);s2.SetRenderColor(cvplot::color::Blue);s1.AddValues(bin1);s2.AddValues(bin2);stringstream ss;ss << "exp_dist lambda:";ss << lambda;ss << " sample:";ss << N;cvplot::View v(ss.str(), { 1200,800 });v.AddSeries(s1);v.AddSeries(s2);cvplot::Figure f;f.SetSize({ 1200,800 });f.SetView(v, 1, 1);f.Show();return 0;
}

结果：