第十五课.马尔科夫链蒙特卡洛方法

在 第十四课中讲述了马尔科夫链与其稳态的性质，本篇讨论基于马尔科夫链蒙特卡洛（MCMC）方法的采样。

M-H采样

Metropolis-Hastings采样原理

我们的目标分布是$p(z)$，同时我们手里有一个便于随时间进行遍历的马尔科夫链，其状态转移矩阵为$Q$。

为了便于在马尔科夫链上随时间进行状态转移，这里的矩阵$Q$设计为：
$$Q_{ij}=P(x^{\ast }|x)=N(\mu =x,\sigma =1).pmf(x^{\ast })$$
上面式子指的是从$t$时刻的某个指定采样点$x$转移到$t+1$时刻的指定采样点$x^{\ast }$的概率等于在均值为$x$，方差为1的整体分布中取得$x^{\ast }$的概率值，本应写作$pdf$，因为正态分布是连续型分布，但是由于是机器采样，实际上已经把目标分布以及正态分布都离散化了，所以$pmf$符合实际操作。

但这样一个（提议矩阵）$Q$显然和目标分布$p(z)$不能满足马尔科夫链的细致平稳条件，即：
$$p(z)Q(z,z^{\ast })\ne p(z^{\ast })Q(z^{\ast },z)$$
这里，为了表达的统一，将状态转移矩阵对应的马尔科夫链的状态$z$转移到状态$z^{\ast }$的概率记作 $Q(z,z^{\ast })$；

为了满足细致平稳条件，引入接受率因子的表达式：
$$\alpha (z,z^{\ast })=min(1,\frac{p(z^{\ast })Q(z^{\ast },z)}{p(z)Q(z,z^{\ast })})$$
原来的不等式两边各自乘上接受率因子，细致平稳条件得到满足：
$$p(z)Q(z,z^{\ast })\alpha (z,z^{\ast })=p(z^{\ast })Q(z^{\ast },z)\alpha (z^{\ast },z)$$
证明如下：

但是，现在的$p(z)$已经不是$Q$对应的马尔科夫链的稳态分布，而是叠加了接受率之后的新转移过程$Q\alpha$。

---

细致平稳条件的意义在于证明一个分布确实是该转移过程对应的稳态分布。

---

M-H采样步骤

假设要采样$N$个样本，初始化第一个时刻的状态$z^{(1)}$，根据马尔科夫链生成$z^{(1)},...,z^{(T)}$；

共$T$轮采样，第$i$轮的采样过程如下：

第一步：从$[0,1]$均匀分布中随机抽取一个数$u$，即$u\sim U(0,1)$；

第二步：利用状态转移矩阵，从上一轮的样本$z^{(i-1)}$依概率抽取出这一步应该转移到的状态$z^{\ast }$，即：$$z^{\ast }\sim Q(z,z^{(i-1)})$$

第三步：计算接受率：$$\alpha =min(1,\frac{p(z^{\ast })Q(z^{\ast },z^{(i-1)})}{p(z^{(i-1)})Q(z^{(i-1)},z^{\ast })})$$
其中，$p(z^{\ast })$可由状态$z^{\ast }$出现的频次统计得到；

第四步：如果$u\le \alpha$，那么我们把$z^{\ast }$作为第$i$轮的采样值，即$z^{(i)}=z^{\ast }$，如果$u\ge \alpha$，抛弃$z^{\ast }$，沿用上一轮的采样值$z^{(i-1)}$作为本轮的采样值，即$z^{(i)}=z^{(i-1)}$

---

其中，第二步的采样体现了转移过程$Q$，第三步和第四步的采样基于$Q$的采样结果，实现了转移过程$Q\alpha$的采样；

---

经过燃烧期后，我们取连续$N$个时间点的状态值$z^{(t)},...,z^{(t+N)}$，作为目标分布的样本近似。

一般而言，我们会选取一个比较适合进行状态转移操作的马尔科夫链，这里描述为：作为状态空间$S$里的任意两个状态，第$i$个状态$x$，第$j$个状态$x^{\ast }$，$Q_{ij}=P(x^{\ast }|x)=N(\mu =x,\sigma =1).pmf(x^{\ast })$；

显然，这样的$Q$是合适的：

• 第一，矩阵的任意一行之和为1，符合概率分布的定义：$$\sum _j{Q}_{ij}=\sum _{x^{\ast }\in S}N(\mu =x,\sigma =1).pmf(x^{\ast })=1$$
• 采样方便，从当前状态$x$，决定下一个状态$x^{\ast }$，就是从正态分布$N(\mu =x,\sigma =1)$中采样一个样本的过程。
• 具有对称性，$Q_{ij}=Q_{ji}$，使得接受率得到化简：$$\alpha =min(1,\frac{p(z^{\ast })Q(z^{\ast },z)}{p(z)Q(z,z^{\ast })})=min(1,\frac{p(z^{\ast })}{p(z)})$$
  从而方便计算。

Gibbs方法

Gibbs核心流程

在M-H的基础上，介绍Gibbs采样；

吉布斯采样是一种针对高维分布的采样方法，假设我们待采样的高维目标分布为$p(z_1,z_2,...,z_m)$，吉布斯采样的目标就是从这个高维分布中采样出一组样本$z^{(1)},z^{(2)},...,z^{(N)}$，使得这组样本服从$m$维分布$p(z)$；

我们约定标记：按照上述假设，从分布中采样得到的样本$z$是一个$m$维的样本，$z_i$为该样本的第$i$维属性值，$z_{-i}$表示除第$i$维外的所有$m-1$个属性值。

Gibbs的核心为一维一维地采样：

以三维随机变量$z$为例：$p(z)=p(z_1,z_2,z_3)$

首先在第一个时刻即时刻0，给定初值：$z_1^{(0)},z_2^{(0)},z_3^{(0)}$

在下一时刻，采样值生成为：$$z_1^{(1)}\sim p(z_1|z_2^{(0)},z_3^{(0)})$$$$z_2^{(1)}\sim p(z_2|z_1^{(1)},z_3^{(0)})$$$$z_3^{(1)}\sim p(z_3|z_1^{(1)},z_2^{(1)})$$

接着往下递推，有：
$$z_1^{(t+1)}\sim p(z_1|z_2^{(t)},z_3^{(t)})$$$$z_2^{(t+1)}\sim p(z_2|z_1^{(t+1)},z_3^{(t)})$$$$z_3^{(t+1)}\sim p(z_3|z_1^{(t+1)},z_2^{(t+1)})$$
其中，$p(z_i|z_{-i})$可以由统计结果得到；

依次循环，采样出$N$个样本：$z^{(1)},z^{(2)},...,z^{(N)}$，去除燃烧期的样本，剩下的就是目标的高维分布。

Gibbs采样的合理性证明

Gibbs其实是M-H采样的特殊情况，首先，我们采样的目标分布为$p(z)$，并把M-H使用的转移矩阵$Q$重新定义为$p$本身：$$Q(z^{\ast },z)=p(z_i|z_{-i}^{\ast })$$；

重新计算接受率：$$\alpha =min(1,\frac{p(z^{\ast })Q(z^{\ast },z)}{p(z)Q(z,z^{\ast })})$$注意到$p(z)=p(z_i,z_{-i})=p(z_i|z_{-i})p(z_{-i})$，所以接受率为：$$\alpha =min(1,\frac{p(z^{\ast })Q(z^{\ast },z)}{p(z)Q(z,z^{\ast })})=min(1,\frac{[p(z_i^{\ast }|z_{-i}^{\ast })p(z_{-i}^{\ast })]p(z_i|z_{-i}^{\ast })}{[p(z_i|z_{-i})p(z_{-i})]p(z_i^{\ast }|z_{-i})})$$
在$z$到$z^{\ast }$的转移中，固定采样值$z$的$-i$维属性，单一采样$z^{\ast }$的第$i$维属性，因此，在一次属性采样前后的$z_{-i}$和$z_{-i}^{\ast }$是等价的，所以将式子中的$z_{-i}^{\ast }$替换为$z_{-i}$，即：

即有Gibbs采样的接受率为1，所以Gibbs采样中不需要像M-H一样生成$u$，并用$u$与接受率进行比较（因为$u$总是小于$\alpha$，所以$z^{\ast }$必然不会被抛弃）。另外，Gibbs是M-H的特例（转移矩阵为目标分布的条件概率，接受率恒为1），所以Gibbs采样的$p$和$Q$满足细致平稳条件，因此按照$Q$的转移过程进行采样，即可以近似得到目标分布$p(z)$。

Gibbs采样实验

令目标分布是一个二维高斯分布$z\sim N(\mu ,\Sigma )$，均值向量为$\mu =[3,3{]}^T$，协方差矩阵为：

现在利用Gibbs采样，获取服从上述分布的样本值。

为了简单演示，我们用已知的条件概率直接生成样本，写成条件分布形式为：

则条件分布$z_1|z_2$和$z_2|z_1$服从：

采样实验如下：

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt


#依照z1|z2的条件高斯分布公式，给定z2的条件情况下采样出z1
def p_z1_given_z2(z2, mu, sigma):
    mu = mu[0] + sigma[1][0] / sigma[0][0] * (z2 - mu[1])
    sigma = sigma[0][0] - sigma[1][0]**2 / sigma[1][1]
    return np.random.normal(mu, sigma)

#依照z2|z1的条件高斯分布公式，给定z1的条件情况下采样出z2
def p_z2_given_z1(z1, mu, sigma):
    mu = mu[1] + sigma[0][1] / sigma[1][1] * (z1 - mu[0])
    sigma = sigma[1][1] - sigma[0][1]**2 / sigma[0][0]
    return np.random.normal(mu, sigma)

#Gibbs采样过程
def gibbs_sampling(mu, sigma, samples_period):
    samples = np.zeros((samples_period, 2))
    
    # np.random.rand()：从0-1均匀分布采样1个数
    z2 = np.random.rand() * 10

    for i in range(samples_period):
        z1 = p_z1_given_z2(z2, mu, sigma)
        z2 = p_z2_given_z1(z1, mu, sigma)
        samples[i, :] = [z1, z2]

    return samples


#目标分布p(z)
mus = np.array([3, 3])
sigmas = np.array([[1, 0.6], [0.6, 1]])

#确定总的采样期和燃烧期
burn_period = int(1e4)
samples_period = int(1e5)
#得到采样样本，舍弃燃烧期样本点
samples = gibbs_sampling(mus, sigmas, samples_period)[burn_period:]
plt.plot(samples[:, 0], samples[:, 1], 'ro', alpha=0.05, markersize=1)
plt.show()