MCMC 浅谈

1. 采样（sampling）是什么

MCMC在采样算法中有着举足轻重的地位，那么什么是采样？采样就是根据某种分布生成样本。举个例子，线性同余发生器就是根据均匀分布生成样本，这就很简单的采样了。

2. 蒙特卡罗算法

假设现在我们有个如下所示的定积分需要进行计算，

\[\theta = \int_{a}^{b}f(x)dx \]

但是$f(x)$的原函数却非常难以获得。这时候采用蒙特卡罗算法进行随机模拟可以得到近似解。假设我们的函数图像如下图所示

我们可以在$[a,b]$之间随机生成一个点$x_0$，然后可以得到$f(x_0)$。假设$[a,b]$间的所有结果都是$f(x_0)$，那么上式的积分结果可以用下式表达

\[\int_a^bf(x)dx \approx (b - a)f(x_0) \]

此时的$f(x_0)$实际上相当于高，$b - a$就相当于底边。用$f(x_0)$表示$[a,b]$区间所有的值显然过于粗糙，那么我们可以生成$n$个样本，用这$n$个样本的均值来计算积分。因为积分是在计算阴影部分的面积，那么积分的结果可以表示成下式

\[\int_a^bf(x)dx \approx \dfrac{(b - a)}{n}\sum_{i=0}^{n-1}f(x_i) \]

此时的$\dfrac{\sum_{i=0}^{n-1}f(x_i)}{n}$可以看成是高。$n$越大，模拟得到的结果就更加贴近真实积分值。

虽然上面的方法可以一定程度上求解出近似的解，但是它隐含了一个假定，即$x$在$[a,b]$之间是均匀分布的，而绝大部分情况，$x$在$[a,b]$之间不是均匀分布的。如果我们用上面的方法，则模拟求出的结果很可能和真实值相差甚远。如果我们可以得到$x$$在[a,b]$的概率分布函数$p(x)$，那么我们的定积分求和可以这样进行：

\[\theta =\int_a^bf(x)dx=\int_a^b\dfrac{f(x)}{p(x)}p(x)dx = \dfrac{1}{n}\sum_{i=0}^{n-1}\dfrac{f(x_i)}{p(x_i)} \]

上式最右边的这个形式就是蒙特卡罗方法的一般形式。当然这里是连续函数形式的蒙特卡罗方法，但是在离散时一样成立。可以看出，最上面我们假设$x$在$[a,b]$之间是均匀分布的时候，$p(x_i)=\dfrac{1}{(b−a)}$，代入我们有概率分布的蒙特卡罗积分的上式，可以得到：

\[\dfrac{1}{n}\sum_{i=0}^{n-1}\dfrac{f(x_i)}{p(x_i)} = \dfrac{1}{n}\sum_{i=0}^{n-1}\dfrac{f(x_i)}{\dfrac{1}{(b−a)}} = \dfrac{(b - a)}{n}\sum_{i=0}^{n-1}f(x_i) \]

3. 接受-拒绝采样

我们常见的概率分布，无论是连续的还是离散的分布，都可以基于Uniform(0,1) 的样本生成。例如正态分布可以通过著名的 Box-Muller 变换得到。其它几个著名的连续分布，包括指数分布、Gamma 分布、t 分布、F 分布、Beta 分布、Dirichlet 分布等等,也都可以通过类似的数学变换得到；离散的分布通过均匀分布更加容易生成。

对于概率分布不是常见的分布，一个可行的办法是采用接受-拒绝采样来得到该分布的样本。既然 $p(x)$ 太复杂在程序中没法直接采样，那么我设定一个程序可采样的分布 $q(x)$，比如高斯分布，然后按照一定的方法拒绝某些样本，以达到接近 $p(x)$ 分布的目的，其中$q(x)$叫做 proposal distribution，$p(x)$叫做target distribution。

接受-拒绝采样

具体采用过程如下，设定一个方便采样的常用概率分布函数 $q(x)$，以及一个常量 $k$，使得 $p(x)$ 总在 $kq(x)$ 的下方，如上图。首先，采样得到$q(x)$的一个样本$z_0$，采样方法如第三节。然后，从均匀分布$(0,kq(z_0))$中采样得到一个值$u$。如果$u$落在了上图中的灰色区域，则拒绝这次抽样，否则接受这个样本$z_0$。重复以上过程得到$n$个接受的样本$z_0,z_1,\dots z_{n−1}$,则最后的蒙特卡罗方法求解结果为：

\[\dfrac{1}{n}\sum_{i=0}^{n-1}\dfrac{f(x_i)}{p(x_i)} \]

整个过程中，我们通过一系列的接受拒绝决策来达到用$q(x)$模拟$p(x)$概率分布的目的。

使用接受-拒绝采样，我们可以解决一些概率分布不是常见的分布的时候，得到其采样集并用蒙特卡罗方法求和的目的。但是接受-拒绝采样也只能部分满足我们的需求，在很多时候我们还是很难得到我们的概率分布的样本集。比如：

对于一些二维分布$p(x,y)$，有时候我们只能得到条件分布$p(x|y)$和$p(y|x)$和,却很难得到二维分布$p(x,y)$一般形式，这时我们无法用接受-拒绝采样得到其样本集。
对于一些高维的复杂非常见分$布p(x_1,x_2,\dots,x_n)$，我们要找到一个合适的$q(x)$和$k$非常困难。

从上面可以看出，要想将蒙特卡罗方法作为一个通用的采样模拟求和的方法，必须解决如何方便得到各种复杂概率分布的对应的采样样本集的问题。下面要介绍的马尔科夫链为解决该问题提供了可行性方案。

4. 马氏链及其平稳分布

（一阶）马尔科夫链定义本身比较简单，它假设某一时刻状态转移的概率只依赖于它的前一个状态。　如果用精确的数学定义来描述，则假设我们的序列状态是$...,X_{t−2},X_{t−1},X_t,X_{t+1},...$，那么我们的在时刻$X_{t+1}$的状态的条件概率仅仅依赖于时刻$X_t$，即：

\[P(X_{t+1} |...X_{t-2}, X_{t-1}, X_{t} ) = P(X_{t+1} | X_{t}) \]

我们先来看马氏链的一个具体的例子。社会学家经常把人按其经济状况分成3类：下层(lower-class)、中层(middle-class)、上层(upper-class)，我们用1,2,3 分别代表这三个阶层。社会学家们发现决定一个人的收入阶层的最重要的因素就是其父母的收入阶层。如果一个人的收入属于下层类别，那么他的孩子属于下层收入的概率是 0.65, 属于中层收入的概率是 0.28, 属于上层收入的概率是 0.07。事实上，从父代到子代，收入阶层的变化的转移概率如下

使用矩阵的表示方式，转移概率矩阵记为

\[P = \begin{bmatrix} 0.65 & 0.28 & 0.07 \\ 0.15 & 0.67 & 0.18 \\ 0.12 & 0.36 & 0.52 \\ \end{bmatrix} \]

假设当前这一代人处在下层、中层、上层的人的比例是概率分布向量 $\pi_0=[\pi_0(1), \pi_0(2), \pi_0(3)]$，那么他们的子女的分布比例将是$\pi_2 = \pi_0 P$,他们的孙子代的分布比例将是 $\pi_2 = \pi_1 P = \pi_0 P^2$,第n代子孙的收入分布比例将是$\pi_n = \pi_{n-1}P = \pi_0P^n$。我们随机生成$\pi_0$，计算各代的结果，代码如下

import numpy as np


def main():
    state = np.random.uniform(size=3)[:, None]
    state = state / np.sum(state, axis=0, keepdims=True)
    
    P = np.array([[0.65, 0.28, 0.07], [0.15, 0.67, 0.18], [0.12, 0.36, 0.52]])

    for i in range(40):
        state = np.dot(P.T, state)
        if i > 20:
            print(state)
            print("*" * 10)

if __name__ == "__main__":
    main()

我们发现给定不同的初始概率分布，最终都收敛到 $\pi=[0.286,0.489,0.225]$这个概率分布。自然的，这个收敛现象并非是我们这个马氏链独有的，而是绝大多数马氏链的共同行为，关于马氏链的收敛我们有如下漂亮的定理：

马氏链定理：如果一个非周期马氏链具有转移概率矩阵$P$，且它的任何两个状态是连通的，那么$\lim_{n->\infty}P_{ij}^n$存在且与$i$无关。记$\lim_{n->\infty}P_{ij}^n = \pi(j)$，我们有

\[\lim_{n\rightarrow\infty}P^n = \begin{bmatrix}\pi(1) & \pi(2) & \cdots & \pi(j) & \cdots \\\pi(1) & \pi(2) & \cdots & \pi(j) & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots\\ \pi(1) & \pi(2) & \cdots & \pi(j) & \cdots \\\cdots & \cdots & \cdots & \cdots & \cdots \end{bmatrix} \]

其中

\[\pi(j) = \sum_{i=0}^\infty \pi(i)P_{ij} $$, $\pi$是方程$\pi = \pi P$的唯一非负解，其中 \]

\pi = [\pi(1), \pi(2), \cdots, \pi(j), \cdots], \sum_{i=1}^\infty \pi(i) = 1

\[$\pi$称为马氏链的平稳分布。代码如下 ```python3 import numpy as np def main(): state = np.random.uniform(size=3)[:, None] state = state / np.sum(state, axis=0, keepdims=True) P = np.array([[0.65, 0.28, 0.07], [0.15, 0.67, 0.18], [0.12, 0.36, 0.52]]) s = np.array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]) for i in range(2, 30): s = np.dot(P,s) if i > 20: print(s) print("*" * 10) print(np.dot(s.T, state)) if __name__ == "__main__": main() ``` 这个马氏链的收敛定理非常重要，**所有的 MCMC(Markov Chain Monte Carlo) 方法都是以这个定理作为理论基础的**。定理的证明相对复杂，一般的随机过程课本中也不给证明，所以我们就不用纠结它的证明了，直接用这个定理的结论就好了。我们对这个定理的内容做一些解释说明： 1. 该定理中马氏链的状态不要求有限，可以是有无穷多个的； 2. 定理中的“非周期“这个概念我们不打算解释了，因为我们遇到的绝大多数马氏链都是非周期的； 3. 两个状态$i,j$是连通并非指$i$ 可以直接一步转移到$j(P_{ij}>0)$,而是指 $i$ 可以通过有限的$n$步转移到达$j(P^n_{ij}>0)$。马氏链的任何两个状态是连通的含义是指存在一个$n$, 使得矩阵$P^n$ 中的任何一个元素的数值都大于零。如果我们得到了某个平稳分布所对应的马尔科夫链状态转移矩阵，我们就很容易采样出这个平稳分布的样本集。假设我们任意初始的概率分布是$\pi_0(x)$, 经过第一轮马尔科夫链状态转移后的概率分布是$\pi_1(x)$，。。。第i轮的概率分布是$\pi_i(x)$。假设经过n轮后马尔科夫链收敛到我们的平稳分布$\pi(x)$，即： \]

\pi_n(x) = \pi_{n+1}(x) = \pi_{n+2}(x) =... = \pi(x)

\[对于每个分布$\pi_i(x)$，我们有：$$\pi_i(x) = \pi_{i-1}(x)P = \pi_{i-2}(x)P^2 = \pi_{0}(x)P^i\]

　
现在我们可以开始采样了，首先，基于初始任意简单概率分布比如高斯分布$\pi_0(x)$采样得到状态值$x_0$，基于条件概率分布$P(x|x_0)$采样状态值$x_1$，一直进行下去，当状态转移进行到一定的次数时，比如到$n$次时，我们认为此时的采样集$(x_n,x_{n+1},x_{n+2},\dots)$即是符合我们的平稳分布的对应样本集。

5. MCMC & Metropolis-Hastings

定理：[细致平稳条件] 如果非周期马氏链的转移矩阵$P$和$\pi(x)$满足

\[\pi(i)P_{ij} = \pi(j)P_{ji} \]

则 $\pi(x)$ 是马氏链的平稳分布。上式被称为细致平稳条件(detailed balance condition)。

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

假设我们已经有一个转移矩阵为$Q$马氏链，$q(i,j)$表示从状态i转移到状态$j$的概率，也可以写为 $q(j|i)$或者$q(i→j)$, 显然，通常情况下

\[p(i) q(i,j) \neq p(j) q(j,i) \]

也就是细致平稳条件不成立，所以 $p(x)$ 不太可能是这个马氏链的平稳分布。我们可否对马氏链做一个改造，使得细致平稳条件成立呢？譬如，我们引入一个 $\alpha(i,j)$, 我们希望

\[p(i) q(i,j)\alpha(i,j) = p(j) q(j,i)\alpha(j,i) $$最简单的，按照对称性，当 $$\alpha(i,j)=p(j)q(j,i)，\alpha(j,i)=p(i)q(i,j)\]

上式便成立了。所以有

\[p(i) \underbrace{q(i,j)\alpha(i,j)}_{Q'(i,j)} = p(j) \underbrace{q(j,i)\alpha(j,i)}_{Q'(j,i)} \]

于是我们把原来具有转移矩阵$Q$的一个很普通的马氏链，改造为了具有转移矩阵$Q′$的马氏链，而 $Q′$恰好满足细致平稳条件，由此马氏链$Q′$的平稳分布就是$p(x)$！

在改造 $Q$的过程中引入的 $\alpha(i,j)$称为接受率，物理意义可以理解为在原来的马氏链上，从状态 $i$ 以$q(i,j)$ 的概率转跳转到状态$j$ 的时候，我们以$\alpha(i,j)$的概率接受这个转移，于是得到新的马氏链$Q′$的转移概率为$q(i,j)\alpha(i,j)$。

假设我们已经有一个转移矩阵$Q$(对应元素为$q(i,j)$), 把以上的过程整理一下，我们就得到了如下的用于采样概率分布$p(x)$的算法。

MCMC算法

上述过程中 $p(x),q(x|y)$ 说的都是离散的情形，事实上即便这两个分布是连续的，以上算法仍然是有效，于是就得到更一般的连续概率分布 $p(x)$的采样算法，而 $q(x|y)$ 就是任意一个连续二元概率分布对应的条件分布。

以上的 MCMC 采样算法已经能很漂亮的工作了，不过它有一个小的问题：马氏链Q在转移的过程中的接受率 $\alpha(i,j)$ 可能偏小，这样采样过程中马氏链容易原地踏步，拒绝大量的跳转，这使得马氏链遍历所有的状态空间要花费太长的时间，收敛到平稳分布$p(x)$的速度太慢。有没有办法提升一些接受率呢?

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

\[p(i)q(i,j)\times 0.1 = p(j)q(j,i) \times 0.2 \]

上式两边扩大5倍，我们改写为

\[p(i)q(i,j) \times 0.5 = p(j)q(j,i) \times 1\]

看，我们提高了接受率，而细致平稳条件并没有打破！这启发我们可以把细致平稳条件式的$\alpha(i,j),\alpha(j,i)$ 同比例放大，使得两数中最大的一个放大到1，这样我们就提高了采样中的跳转接受率。所以我们可以取$$\alpha(i,j) = \min\left{\frac{p(j)q(j,i)}{p(i)q(i,j)},1\right}$$

于是，经过对上述MCMC 采样算法中接受率的微小改造，我们就得到了如下教科书中最常见的 Metropolis-Hastings 算法。
Metropolis-Hastings算法

　很多时候，我们选择的马尔科夫链状态转移矩阵$Q$如果是对称的，即满足$Q(i,j)=Q(j,i)$,这时我们的接受率可以进一步简化为：$$\alpha(i,j) = min{ \frac{\pi(j)}{\pi(i)},1}$$
　
为了更容易理解，这里给出一个M-H采样的实例。我们的目标平稳分布是一个均值3，标准差2的正态分布，而选择的马尔可夫链状态转移矩阵Q(i,j)的条件转移概率是以i为均值,方差1的正态分布在位置j的值。

import random
import math
import matplotlib
from scipy.stats import norm
matplotlib.use('Qt5Agg')
import matplotlib.pyplot as plt



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


def main(sigma=1, num_bins=50, T=50000, t = 45000):
    pi = [0]
    for _ in range(T):
        pi_star = norm.rvs(loc=pi[-1], scale=sigma, size=1, random_state=None)[0]
        alpha = min(1, norm_dist_prob(pi_star)/norm_dist_prob(pi[-1]))
        pi.append(pi_star if random.uniform(0, 1) < alpha else pi[-1])
    
    plt.scatter(pi[-t:], norm.pdf(pi[-t:], loc=3, scale=2))
    plt.hist(pi[-t:], num_bins, normed=1, facecolor='red', alpha=0.7)
    plt.show()

if __name__ == "__main__":
    main()

M-H采样完整解决了使用蒙特卡罗方法需要的任意概率分布样本集的问题，因此在实际生产环境得到了广泛的应用。但是在大数据时代，M-H采样面临着两大难题：

我们的数据特征非常的多，M-H采样由于接受率计算式$\dfrac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)}$的存在，在高维时需要的计算时间非常的可观，算法效率很低。同时$\alpha(i,j)$一般小于1，有时候辛苦计算出来却被拒绝了。能不能做到不拒绝转移呢？
由于特征维度大，很多时候我们甚至很难求出目标的各特征维度联合分布，但是可以方便求出各个特征之间的条件概率分布。这时候我们能不能只有各维度之间条件概率分布的情况下方便的采样呢？

接下来要讲的Gibbs采样很好地解决了上面两个问题。

6. Gibbs 采样

对于高维的情形，由于接受率 $\alpha$的存在(通常 $\alpha <1$), 以上 Metropolis-Hastings 算法的效率不够高。能否找到一个转移矩阵$Q$使得接受率 $\alpha=1$ 呢？我们先看看二维的情形，假设有一个概率分布 $p(x,y)$, 考察x坐标相同的两个点$A(x_1,y_1),B(x_1,y_2)$，我们发现

\[\begin{align*} p(x_1,y_1)p(y_2|x_1) = p(x_1)p(y_1|x_1)p(y_2|x_1) \\ p(x_1,y_2)p(y_1|x_1) = p(x_1)p(y_2|x_1)p(y_1|x_1) \end{align*} \]

所以得到

\[p(x_1,y_1)p(y_2|x_1) = p(x_1,y_2)p(y_1|x_1) \]

即

\[p(A)p(y_2|x_1) = p(B)p(y_1|x_1) \]

基于以上等式，我们发现，在 $x=x_1$ 这条平行于 y轴的直线上，如果使用条件分布 $p(y|x_1)$做为任何两个点之间的转移概率，那么任何两个点之间的转移满足细致平稳条件。同样的，如果我们在 $y=y_1$ 这条直线上任意取两个点 $A(x_1,y_1),C(x_2,y_1)$,也有如下等式

\[p(A)p(x_2|y_1) = p(C)p(x_1|y_1). \]

平面上马氏链转移矩阵的构造

于是我们可以如下构造平面上任意两点之间的转移概率矩阵$Q$

\[\begin{align*} Q(A\rightarrow B) & = p(y_B|x_1) & \text{if} \quad x_A=x_B=x_1 & \\ Q(A\rightarrow C) & = p(x_C|y_1) & \text{if} \quad y_A=y_C=y_1 & \\ Q(A\rightarrow D) & = 0 & \text{otherwise} & \end{align*} \]

有了如上的转移矩阵$Q$, 我们很容易验证对平面上任意两点 $X,Y$，满足细致平稳条件

\[p(X)Q(X\rightarrow Y) = p(Y) Q(Y\rightarrow X) \]

下面给出吉布斯采样算法的伪代码

以上算法收敛后，得到的就是概率分布$p(x_1,x_2,⋯,x_n)$的样本，当然这些样本并不独立，但是我们此处要求的是采样得到的样本符合给定的概率分布，并不要求独立。同样的，在以上算法中，坐标轴轮换采样不是必须的，可以在坐标轴轮换中引入随机性，这时候转移矩阵 $Q$ 中任何两个点的转移概率中就会包含坐标轴选择的概率，而在通常的 Gibbs Sampling 算法中，坐标轴轮换是一个确定性的过程，也就是在给定时刻$t$，在一根固定的坐标轴上转移的概率是1。

假设我们要采样的是一个二维正态分布$Norm(\mu,\Sigma)$,其中：

\[\begin{split} \mu &= (\mu_1,\mu_2) &= (5,-1) \\ \Sigma &= \left( \begin{array}{ccc} \sigma_1^2&\rho\sigma_1\sigma_2 \\ \rho\sigma_1\sigma_2 &\sigma_2^2 \end{array} \right) &= \left( \begin{array}{ccc} 1&1 \\ 1&4 \end{array} \right) \end{split} \]

而采样过程中的需要的状态转移条件分布为：

\[\begin{split} P(x_1|x_2) &= Norm\left ( \mu _1+\rho \sigma_1/\sigma_2 \left ( x _2-\mu _2 \right ), (1-\rho ^2)\sigma_1^2 \right )\\ P(x_2|x_1) &= Norm\left ( \mu _2+\rho \sigma_2/\sigma_1 \left ( x _1-\mu _1 \right ), (1-\rho ^2)\sigma_2^2 \right ) \end{split} \]

下面给出代码

from mpl_toolkits.mplot3d import Axes3D
from scipy.stats import multivariate_normal
import random
import math
import matplotlib
matplotlib.use('Qt5Agg')
import matplotlib.pyplot as plt


def pxy(xy, m1, m2, s1, s2, rho=0.5):
    return random.normalvariate(m2 + rho * s2 / s1 * (xy - m1), math.sqrt(1 - rho**2)*s2)


def Gibbs(N=5000, K=20, m1=5, m2=-1, s1=1, s2=2, num_bins=50):
    x, y, result = 0, m2, []
    samplesource = multivariate_normal(mean=[5, -1], cov=[[1, 0.5], [0.5, 2]])
    for i in range(N):
        for j in range(K):
            x = pxy(y, m2, m1, s2, s1)
            y = pxy(x, m1, m2, s1, s2)
            z = samplesource.pdf([x, y])
            result.append((x, y, z))
    X, Y, Z = zip(*result)
    plt.hist(X, num_bins, density=1, facecolor='green', alpha=0.5)
    plt.hist(Y, num_bins, density=1, facecolor='red', alpha=0.5)
    plt.title('Histogram')

    fig = plt.figure()
    ax = Axes3D(fig, rect=[0, 0, 1, 1], elev=30, azim=20)
    ax.scatter(X, Y, Z, marker='o')
    plt.show()


if __name__ == "__main__":
    Gibbs()

输出的两个特征各自的分布如下：