MCMC采样及M-H采样

前言

前面的文章已经说到给定一个概率平稳分布π, 只要能够其马尔科夫链状态转移矩阵P，我们就可以找到一种通用的概率分布采样方法，进而用于蒙特卡罗模拟。现在就学习下解决这个问题的办法：MCMC采样和它的易用版M-H采样。

如何能做到这一点呢？我们主要使用如下的定理。

马尔科夫链的细致平稳条件

细致平稳条件: 如果非周期马尔科夫链的状态转移矩阵$P$和概率分布$π(x)$对于所有的$i,j$满足：

$$\pi(i)P(i,j) = \pi(j)P(j,i)$$

则称概率分布$π(x)$是状态转移矩阵$P$的平稳分布。上式被称为细致平稳条件(detailed balance condition)。

其实这个定理是显而易见的，因为细致平稳条件的物理含义就是对于任何两个状态$i,j$, 从$i$转移出去到$j$而丢失的概率质量，恰好会被从$j$转移回$i$的概率质量补充回来，所以状态上的概率质量$π(i$)是稳定的，从而$π(x)$是马尔科夫链的平稳分布。

证明很简单,由细致平稳条件可知:

$$\sum\limits_{i=1}^{\infty}\pi(i)P(i,j) = \sum\limits_{i=1}^{\infty} \pi(j)P(j,i) = \pi(j)\sum\limits_{i=1}^{\infty} P(j,i) = \pi(j)$$

将上式用矩阵表示即为：

$$\pi P = \pi$$

即满足马尔可夫链的收敛性质。也就是说，只要我们找到了可以使概率分布$π(x)$满足细致平稳分布的矩阵$P$即可。这给了我们寻找从平稳分布$π$, 找到对应的马尔科夫链状态转移矩阵$P$的新思路。

不过不幸的是，仅仅从细致平稳条件还是很难找到合适的矩阵P。比如我们的目标平稳分布是$π(x)$,随机找一个马尔科夫链状态转移矩阵$Q$,它是很难满足细致平稳条件的，即:

$$\pi(i)Q(i,j) \neq \pi(j)Q(j,i)$$

那么如何使这个等式满足呢？

MCMC采样

由于一般情况下，目标平稳分布$π(x)$和某一个马尔科夫链状态转移矩阵$Q$不满足细致平稳条件，即:

$$\pi(i)Q(i,j) \neq \pi(j)Q(j,i)$$

我们可以对上式做一个改造，使细致平稳条件成立。方法是引入一个$α(i,j)$,使上式可以取等号，即：

$$\pi(i)Q(i,j)\alpha(i,j) = \pi(j)Q(j,i)\alpha(j,i)$$

此时的问题转化为什么样的$α(i,j)$可以使等式成立, 其实很简单，只要满足下两式即可：

$$\alpha(i,j) = \pi(j)Q(j,i) \\ \alpha(j,i) = \pi(i)Q(i,j)$$

这样，我们就得到了我们的分布$π(x)$对应的马尔科夫链状态转移矩阵$P$，满足：

$$P(i,j) = Q(i,j)\alpha(i,j)$$

也就是说，我们的目标矩阵$P$可以通过任意一个马尔科夫链状态转移矩阵$Q$乘以$α(i,j)$得到。$α(i,j)$我们有一般称之为接受率。取值在[0,1]之间，可以理解为一个概率值。即目标矩阵$P$可以通过任意一个马尔科夫链状态转移矩阵$Q$以一定的接受率获得。这个很像我们在蒙特卡罗方法讲到的接受-拒绝采样，那里是以一个常用分布通过一定的接受-拒绝概率得到一个非常见分布，这里是以一个常见的马尔科夫链状态转移矩阵$Q$通过一定的接受-拒绝概率得到目标转移矩阵$P$,两者的解决问题思路是类似的。

MCMC的采样过程如下

1. 输入我们任意选定的马尔科夫链状态转移矩阵$Q$，平稳分布$π(x)$，设定状态转移次数阈值$n_1$，需要的样本个数$n_2$
2. 从任意简单概率分布采样得到初始状态值$x_0$
3. $for\; t=0 \; to \; n_1+n_2−1:$
   1. 从条件概率分布$Q(x|x_t)$中采样得到样本$x_∗$
   2. 从均匀分布采样$u∼uniform[0,1]$
   3. 如果$u < \alpha(x_t,x_{*}) = \pi(x_{*})Q(x_{*},x_t)$, 则接受转移$x_t \to x_{*}$，即$x_{t+1}= x_{*}$
   4. 否则不接受转移，$t=max(t−1,0)$
4. 样本集$(x_{n_1}, x_{n_1+1},..., x_{n_1+n_2-1})$即为我们需要的平稳分布对应的样本集。

以上的MCMC采样算法已经能很漂亮的工作了，不过它有一个小的问题：马尔科夫链$Q$在转移的过程中的接受率$α(i,j)$可能偏小，这样采样过程中容易原地踏步，拒绝大量的跳转，使得马尔科夫链收敛到平稳分布$p(x)$的速度太慢。有可能我们采样了上百万次马尔可夫链还没有收敛，也就是上面这个$n_1$要非常非常的大，这让人难以接受，怎么办呢？

看下节.

M-H采样

M-H采样是Metropolis-Hastings采样的简称，这个算法首先由Metropolis提出，被Hastings改进，因此被称之为Metropolis-Hastings采样或M-H采样.

M-H采样解决了我们上一节MCMC采样接受率过低的问题。

假设$α(i,j)=0.1,α(j,i)=0.2$此时满足细致平稳条件，于是:

$$p(i)q(i,j) \cdot 0.1=p(j)q(j,i) \cdot 0.2$$

上式两边同时扩大5倍，等式变为:

$$p(i)q(i,j) \cdot 0.5=p(j)q(j,i) \cdot 1$$

我们提高了接受率，而细致平稳条件并没有打破。这启发我们可以把细致平稳条件中$α(i,j),α(j,i)$等比例放大，使得两数中较大的一个放大到1，如此提高了采样中的跳转接受率。 这样我们的接受率可以做如下改进，即：

$$\alpha(i,j) = min\{ \frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)},1\}$$

通过这个微小的改造，我们就得到了可以在实际应用中使用的M-H采样算法过程如下：

1. 输入我们任意选定的马尔科夫链状态转移矩阵$Q$，平稳分布$π(x)$，设定状态转移次数阈值$n_1$，需要的样本个数$n_2$
2. 从任意简单概率分布采样得到初始状态值$x_0$
3. $for\; t=0 \;to\; n_1+n_2−1:$
   1. 从条件概率分布$Q(x|x_t)$中采样得到样本$x_∗$
   2. 从均匀分布采样$u∼uniform[0,1]$
   3. 如果$u < \alpha(x_t,x_{*}) = min\{ \frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)},1\}$, 则接受转移$x_t \to x_{*}$，即$x_{t+1}=x_∗$
   4. 否则不接受转移，$t=max(t-1,0)$
4. 样本集$(x_{n_1}, x_{n_1+1},..., x_{n_1+n_2-1})$即为我们需要的平稳分布对应的样本集。

很多时候，我们选择的马尔科夫链状态转移矩阵$Q$如果是对称的，即满足$Q(i,j)=Q(j,i)$,这时我们的接受率可以进一步简化为：

$$\alpha(i,j) = min\{ \frac{\pi(j)}{\pi(i)},1\}$$

Note注意点

• 具体的关于上述两个采样的方法的转移矩阵$Q$的选择，如果是离散的分布，那么只要这个矩阵满足每行概率和为1即可，没有特殊要求。如果是连续分布，则一般会选择正态分布作为$Q$，因为正态分布的联合分布，条件分布，边缘分布都是正态分布，所以后面计算条件概率比较的容易。对采样的好坏一般影响都不大，只能说不同的$Q$对算法的收敛速度会快有慢，只要收敛了，最后的采样结果都是可以的。
• 这里$α(i,j)$是某一个位置的$(i,j)$的接受率，所以α实际上也是一个矩阵，如果按矩阵的写法，就是:$P=Qα^T$, 不要理解为所有的位置都乘以的是同一个α(i,j),这是不对的。

M-H采样实例

我们的目标平稳分布是一个均值3，标准差2的正态分布，而选择的马尔可夫链状态转移矩阵$Q(i,j)$的条件转移概率是以$i$为均值,方差$1$的正态分布在位置$j$的值。这个例子仅仅用来加深对M-H采样过程的理解。毕竟一个普通的一维正态分布用不着去用M-H采样来获得样本。

import random
import math
from scipy.stats import norm
import matplotlib.pyplot as plt
%matplotlib inline

def norm_dist_prob(theta):
    y = norm.pdf(theta, loc=3, scale=2)
    return y

T = 100000
pi = [0 for i in range(T)]
sigma = 1
t = 0
while t < T-1:
    t = t + 1
    pi_star = norm.rvs(loc=pi[t - 1], scale=sigma, size=1, random_state=None)
    alpha = min(1, (norm_dist_prob(pi_star[0]) / norm_dist_prob(pi[t - 1])))

    u = random.uniform(0, 1)
    if u < alpha:
        pi[t] = pi_star[0]
    else:
        pi[t] = pi[t - 1]

plt.scatter(pi, norm.pdf(pi, loc=3, scale=2))
num_bins = 1000
plt.hist(pi, num_bins, density=1, facecolor='red', alpha=0.7)
plt.show()

结果如下：

M-H采样总结

M-H采样完整解决了使用蒙特卡罗方法需要的任意概率分布样本集的问题，因此在实际生产环境得到了广泛的应用。

但是在大数据时代，M-H采样面临着两大难题：

1. 我们的数据特征非常的多，M-H采样由于接受率计算式$\frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)}$的存在，在高维时需要的计算时间非常的可观，算法效率很低。同时$α(i,j)$一般小于1，有时候辛苦计算出来却被拒绝了。能不能做到不拒绝转移呢？

2. 由于特征维度大，很多时候我们甚至很难求出目标的各特征维度联合分布，但是可以方便求出各个特征之间的条件概率分布。这时候我们能不能只有各维度之间条件概率分布的情况下方便的采样呢？

Gibbs采样解决了上面两个问题，因此在大数据时代，MCMC采样基本是Gibbs采样的天下，下一篇我们就来学习Gibbs采样。