Metropolis-Hasting抽样算法

1. 随机模拟的基本思想

假设我们有一个矩形区域 $R$ ，面积为 $S_0$ 。在此区域中，有一个不规则区域 $M$ ，其面积 $S$ 待求。

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

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

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

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

2. 拒绝抽样

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

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

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

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

拒绝抽样的方法是：

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

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

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

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

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

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

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

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

3. Metropolis-Hastings抽样

3.1 引入思想

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

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

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

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

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

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

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

定理——细致平稳条件：

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

证明：
$∑i=1∞π(i)Pij=∑i=1∞π(j)Pji=π(j)∑i=1∞Pji=π(j),∀jsum_{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$
这说明：
$π (x) P = π (x)$
因此 $π (x)$ 是平稳分布。

物理意义：

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

3.3 M-H算法实现

假设我们已经有了一个转移矩阵 $Q$ 。一般情况下，该转移矩阵不满足细致平稳条件：
$P(i)Qij≠P(j)QjiP(i)Q_{ij} not= P(j)Q_{ji}$
我们引入接受率 $α (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)$ 是状态 $i$ 出现的概率，是假设马氏链满足复杂分布时输出状态的概率。
$Q$ 是初始转移概率矩阵，满足任意一个给定的简单分布，便于抽样即可。

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

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

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

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

图解转移过程：

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

算法描述如下：

初始化马氏链状态 $X_0=x_0$
从 $t = 0$ 开始，重复以下步骤：
- 假设该 $t$ 时刻状态为 $X_t=x_t$ ；
- 对 $Q_{x_t,x}$ 抽样，随机得到一个状态 $x_{next}$ ；
- 在 $[0, 1]$ 上均匀抽样，得到一个 $u$ ；
- 若 $u<α(xt,xnext)=P(xnext)Qxnext,xtu<alpha(x_t,x_{next})=P(x_{next})Q_{x_{next},x_t}$ ，则接受转移 $X_{t+1}=x_{next}$ ；
- 否则不转移，即 $X_{t+1}=x_t$ 。

3.4 算法升级

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

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

解释：

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

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

3.5 仿真实验

仿真目标：

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

简化：

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

代码实现：

#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;
}

结果：

Metropolis-Hasting抽样算法

目录

Metropolis-Hasting抽样算法

1. 随机模拟的基本思想

2. 拒绝抽样

3. Metropolis-Hastings抽样

3.1 引入思想

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

3.3 M-H算法实现

3.4 算法升级

3.5 仿真实验

Published by

风君子