Metropolis Hasting 算法

References：

• 马尔可夫链蒙特卡罗算法（MCMC）
• 走进贝叶斯统计（五）—— Metropolis-Hasting 算法

【引入】

假设参数 $\theta$ 的先验分布 $P(\theta)$ 为某定义域 $\theta\in[0,1]$，均值 $\mu=\frac{1}{2}$，正则化系数为 $\alpha$ 的截尾正态分布，后验分布 $P(\theta|X)$ 具备如下形式：

$$p(\theta|x) = \frac{\theta^{x}(1-\theta)^{1-x}\cdot \frac{\alpha}{\sqrt{2\pi}}e^{-\frac{(\theta-\frac{1}{2})^2}{2\sigma^2}}}{\int_{0}^{1} \theta^{x} (1-\theta)^{1-x} \cdot \frac{\alpha}{\sqrt{2\pi}}e^{-\frac{(\theta-\frac{1}{2})^2}{2\sigma^2}} d\theta}$$

在使用蒙特卡罗法对如上例这样的后验分布 $P(\theta|X)$ 进行抽样时，一个难点是如何处理分母的积分：

$$c=\int_{0}^{1} \theta^{x} (1-\theta)^{1-x} \cdot \frac{\alpha}{\sqrt{2\pi}}e^{-\frac{(\theta-\frac{1}{2})^2}{2\sigma^2}} d\theta$$

对于这个问题，一个直观的解决方法是：能否通过某种操作，来消除掉这个积分 $c$ ？

记后验分布为：

$$P(\theta|X) = \frac{P(X|\theta)P(\theta)}{c}$$

此时，通过计算 $\theta$ 的某两个具体的值 $\theta_a$ 与 $\theta_b$ 的相对概率密度，来从式子中抵消掉 $c$，即：

$$\gamma = \frac{p(\theta_a|x)}{p(\theta_b|x)} = \frac{p(\theta_a|x)p(\theta_a)/c}{p(\theta_b|x)p(\theta_b)/c} = \frac{p(\theta_a|x)p(\theta_a)}{p(\theta_b|x)p(\theta_b)}$$

可以发现，在这样的场景下，成功避免了对 $c$ 的计算，但又有了新的问题：如何利用相对概率密度进行抽样？

在大样本的情况下，$\gamma$ 实际上相当于样本中 $\theta_a$ 的个数与 $\theta_b$ 的个数的比值，即：

$$\gamma = \frac{p(\theta_a|x)}{p(\theta_b|x)} = \frac{\#(\theta_a)}{\#(\theta_b)}$$

这时候可以想到一个理想情况下的抽样方案

若有一组样本 $\{\theta_1,\theta_2,\cdots,\theta_n\}$，恰好这其中有 $a$ 个等于 $\theta_a$，有 $b$ 个等于 $\theta_b$，那么当随机得到一个数 $\theta_{new}$，发现 $\theta_{new}=\theta_a$ 时，就能做出如下决策：

• 当发现 $\theta_a$ 的个数过多，导致 $\frac{a}{b}>\gamma$ 时，就倾向认为这个数 $\theta_{new}$ 不取为好
• 当发现 $\theta_a$ 的个数过少，导致 $\frac{a}{b}<\gamma$ 时，就倾向认为将这个数 $\theta_{new}$ 放入样本中，构造 $\{\theta_1,\theta_2,\cdots,\theta_n,\theta_{new}\}$，因为这样能使 $\frac{a}{b}$ 更加逼近 $\gamma$

当然，上述的这种思路十分理想化，因为在实际应用中，已有的样本中 $\theta$ 的取值繁多，在连续分布中甚至不太可能产生相同取值的样本

但这个思路提供了一个良好的解决方案，能够帮助我们在复杂的后验分布中进行抽样，这种思路即 MH 算法（Metropolis-Hasting）算法的雏形

此外，从上述的这种理想的从抽样方法中可以看到，抽样过程是基于现有样本的，也就是说抽样是基于条件概率 $P(\theta|\theta_1,\theta_2,\cdots,\theta_n,X)$ 的，即在这种情况下，取得的样本不独立

【基本思路】

通过上述抽样思路得到的样本并不是相互独立的，因此在 MH 算法中，常将样本想象成一个有序的马尔可夫链，然后在这个序列中，通过当前样本中的最后一个样本来对下一个样本进行抽样

也就是说，在 MH 算法中，对于目标分布 $P(\theta)$，通过采样后得到了样本序列 $\theta_1,\cdots,\theta_n$（$\theta_i$ 表示在第 $i$ 次采样中得到的样本），然后使用已采的样本 $\theta_n$ 来采样 $\theta_{n+1}$

此时，有一个值得注意的点：若刚开始没有进行采样，对于第一个样本 $\theta_1$ 的确定，通常是任意地选取 $P(\theta)$ 定义域内的一个值

综上，利用 MH 算法进行抽样，只需循环以下两个操作：

1. 基于某种分布 $P(\theta)$，利用 $\theta_{last}$ 抽样得到 $\theta_{new}$
2. 通过 $\gamma$ 来决定新样本 $\theta_{n+1}$ 的值

在这里，存在两个问题，即：

1. 如何构建具体的马尔可夫链，即要如何定义转移概率矩阵或转移核函数，以保证 MCMC 的条件成立？
2. 应当基于怎样的一个分布，并结合 $\theta_{last}$ 来得到 $\theta_{new}$ ？

【马尔可夫链的构造与提议分布】

转移核

假设要抽样的概率分布为 $p(x)$，MH 算法采用转移核为 $p(x,x’)$ 的马尔可夫链，即：

$$p(x,x')=q(x,x')\alpha(x,x')$$

其中，$q(x,x’)$ 为提议分布（Proposal Distribution），$\alpha(x,x’)$ 为接受分布（Acceptance Distribution）

提议分布 $q(x,x’)$ 为一个不可约的马尔可夫链的转移核，其概率值恒不为零，同时是一个容易抽样的分布

接受分布 $\alpha(x,x’)$ 为：

$$\alpha(x,x') = \min \{ 1, \frac{p(x')q(x',x)}{p(x)q(x,x')} \}$$

因此，转移核 $p(x,x’)$ 可写为：

$$p(x,x') = \left\{\begin{array}{aa} q(x,x'), & p(x')q(x',x)\geq p(x)q(x,x') \\ q(x',x)\frac{p(x')}{p(x)}, & p(x')q(x',x)< p(x)q(x,x') \end{array}\right.$$

关于马尔可夫链的转移核，见：马尔可夫链的转移概率

提议分布

提议分布 $q(x,x’)$ 有多种可能的形式，这里只介绍两种 MH 算法中常用的形式

Metropolis 选择

Metropolis 选择是 MH 算法最初采用的提议分布，由于其是对称的，也被称为对称提议分布（Symmetric Proposal Distribution），即

$$q(x,x')=q(x',x)$$

此时，接受分布 $\alpha(x,x’)$ 可化简为：

$$\alpha(x,x') = \min \{ 1, \frac{p(x')}{p(x)} \}$$

对称提议分布的一大特点是：当 $x$ 与 $x’$ 接近时，$q(x,x’)$ 的概率高，状态转移在附近点的可能性更大，否则概率低，状态转移在附近点的可能性更小

此外，在 MH 算法中，对称提议分布的一个特例是取条件概率分布，即 $q(x,x’)=p(x’|x)$，此时有：

$$p(x'|x) = p(x|x')$$

通常来说，会选取会选取均匀分布 $U(x-\delta,x+\delta)$ 或正态分布 $N(x,\delta^2)$ 来作为对称提议分布，其中 $\delta$ 是关系到抽样效率的超参数

独立抽样

假设提议分布 $q(x,x’)$ 与当前状态 $x$ 无关，即：

$$q(x,x')=q(x')$$

提议分布的计算按照 $q(x’)$ 独立抽样进行

此时，接受分布 $\alpha(x,x’)$ 可以写为：

$$\alpha(x,x') = \min \{ 1, \frac{p(x')q(x)}{p(x)q(x')} \}$$

令 $w(x’)=\frac{p(x’)}{q(x’)},w(x)=\frac{p(x)}{q(x)}$，则有：

$$\alpha(x,x') = \min \{ 1, \frac{w(x')}{w(x)} \}$$

独立抽样实现简单，但收敛速度较慢，通常选择接近目标分布 $p(x)$ 的分布作为提议分布 $q(x)$

【算法流程】

综上所述，MH 的算法流程如下

输入：抽样的目标分布的密度函数 $p(x)$，函数 $f(x)$

输出：$p(x)$ 的随机样本 $x_{m+1},x_{m+2},\cdots,x_{n}$，函数样本均值 $\hat{f}_{mn}$

参数：收敛步数 $m$，迭代步数 $n$

算法步骤：

Step 1：算法初始化，随机取得或人工定义第一个样本 $x_0$

Step 2：对 $i=1,2,\cdots,n$ 循环执行

1）令 $x=x_{i-1}$ 从提议分布 $q(x,x’)$ 中抽样，抽取候选状态 $x’=x_{i-1}’$

2）计算接受概率

$$\alpha(x,x') = \min \{ 1,\frac{p(x')q(x',x)}{p(x)q(x,x')} \}$$

3）从均匀分布 $U(0,1)$ 中随机抽取一随机数 $u$，若有：

$$u\leq \alpha(x,x')$$

则令 $x_i=x_{i-1}’$，否则，令 $\theta_i=\theta_{i-1}$

Step 3：从样本集合 $\{x_0,x_{1},x_{2},\cdots,x_{n}\}$ 中得到样本集合 $\{x_{m+1},x_{m+2},\cdots,x_{n}\}$，并计算函数样本均值 $\hat{f}_{mn}$

$$\hat{f}_{mn} = \frac{1}{n-m}\sum_{i=m+1}^n f(x_i)$$

【示例】

如下图所示，需要使用 MH 算法对后验分布 $N(10.04,0.44)$ 进行采样，从直方图可以看到，通过 MH 算法选取的样本很好地表示了目标分布的真实形态

在这个例子中，选取了一个非常离谱的初值 $\theta_1=0$，但是从图中可以看到，即使是在这样的情况下，仍然能通过几百步的迭代使采样恢复到正常区间，因此不必过多的担心初始值 $\theta_1$ 的选取问题

此外，在实际的应用中，往往会将前面部分的抽样结果给截断，以保证样本能够很好表示目标分布