Metropolis-Hasting抽样算法

Metropolis-Hasting抽样算法

1. 随机模拟的基本思想

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

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

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

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

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

2. 拒绝抽样

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

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

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

我们先构造比较简单的、容易抽样的高斯分布$q(z)$，乘一个常数$k$，使之满足：
$$kq(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)$采样。

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

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

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

3. Metropolis-Hastings抽样

3.1 引入思想

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

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

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

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

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