13 MCMC——马尔可夫链蒙特卡洛

13.1 MCMC的意义

我们都知道Monte Carlo是一种常用的估计方法，我们经常用这个方法对期望进行随机估计。而Monte Carlo中的一大问题就是采样，要从一个复杂分布中采样是困难的，我们只能通过一些方法，将复杂的分布表示出来，然后进行采样。

我们从三个方面介绍一下MCMC的意义：

1. 采样的动机：
   • 采样本身：我们经常需要从各种概率分布中取得样本。（例如学习一张图片后生成新的图片）
   • 求和或求积分（期望）：$E_{p(x)}[f(x)]\approx \frac{1}{N}{\sum }_{i=1}^{N}f(x^{(i)})$
2. 什么是好的样本：
   • 样本趋向（集中）于高概率区域——能够凸显分布特征
   • 样本之间相互独立
3. 为什么会采样困难：
   • 无法归一化：Partation funtion is intractable——求不出概率（归一化因子是积分形式，维度过高会无法求解）
   • 维度灾难：高维的稀疏性使得直接判定的方法无法操作

13.2 简单采样方法介绍

首先简单介绍一下Monte Carlo方法：Monte Carlo Method是基于采样的随即近似方法。假如存在$N$​个样本$Z^{(1)},\dots ,Z^{(N)}$​，我们就可以通过代入样本求均值求期望：
$$P(Z|X)\to E_{Z|X}[f(Z)]={\int }_{Z}f(Z)\cdot P(Z|X)d_Z\approx \frac{1}{N}\sum _{i=1}^{N}f(Z^{(i)})$$
同时给出一个最基本的条件：我们认为样本$u\backsim U(0,1)$的均匀分布是很好取的。

13.2.1 概率分布采样

条件：我们已知目标分布$p(x)$的pdf，而且可以通过积分求出其cdf。

方法：因为均匀分布$u\backsim U(0,1)$很好取，所以我们随机取$u$作为cdf的值，通过反函数求出样本。

缺点：大多数分布的积分、反函数求不出来，无法实现

13.2.2 Rejection Sampling——拒绝采样

条件：我们已知目标分布$p(x)$的pdf，用$q(x)$表示一个可以通过概率分布采样的简单的pdf，令$\forall Z\to M\cdot q(Z)\ge p(Z)$（带$M$是因为$q$与$p$的面积相同，不带系数无法实现）

方法：

• 假定在样本点$Z^{(i)}$时，$p(Z^{(i)})$的值为$l_1$，$M\cdot q(Z^{(i)})$的值为$l_2$，如图所示
  

• 定义一个变量，接受率$\alpha =\frac{l_1}{l_2}$，取值范围为$\alpha \in [0,1]$。

• 在对$q(Z)$采样的时候，随机在均匀分布中采样：$u\backsim U(0,1)$。若$u\le \alpha$则保留$Z^{(i)}$，反之则删除。

缺点：若$p(Z^{(i)})$和$M\cdot q(Z^{(i)})$的值相差过大，会增加运算量（可能取好多次样都保留不了几个）

13.2.3 Importance Sampling——重要性采样

原理：对取样分布进行变换：
$$\begin{array}{rl}{E}_{p(Z)}[f(Z)]& ={\int }_{Z}f(Z)p(Z){\mathrm{d}}_Z\\ & ={\int }_{Z}f(Z)\frac{p(Z)}{q(Z)}\cdot q(Z){\mathrm{d}}_Z\\ & =E_{q(Z)}\left[f(Z)\frac{p(Z)}{q(Z)}\right]\\ & \approx \frac{1}{N}\sum _{i=1}^N\left[f(Z^{(i)})\frac{p(Z^{(i)})}{q(Z^{(i)})}\right]\end{array}$$
我们将$q(Z)$​定义为一个简单分布。通过对$q(Z)$​进行采样，可以求出$f(Z^{(i)})\frac{p(Z^{(i)})}{q(Z^{(i)})}$​，我们将$\frac{p(Z^{(i)})}{q(Z^{(i)})}$​称为权重，乘上倒数我们便可以得到结果：
$$\frac{1}{N}\sum _{i=1}^{N}f(Z^{(i)})（对p取样）\to \frac{1}{N}\sum _{i=1}^N\left[f(Z^{(i)})\frac{p(Z^{(i)})}{q(Z^{(i)})}\right]（对q取样）$$
缺点：以来$p$、$q$的相似度，若相差过大，则样本位置分布的差别会很大。

推广：Sampling Importance Resampling——重采样：

• 先进行一次Importance Resampling得到$Z^{(i)}$
• 通过将权重$\frac{p(Z^{(i)})}{q(Z^{(i)})}$设置为概率重新进行采样——可以将更多的点聚集在高概率区域

13.3 Markov Chain知识补充

13.3.1 Markov Chain定义

Markov Chain定义：通过 { x t } {\lbrace x_t \rbrace} {xt​}表示节点，$P=[P_{ij}]$表示转移矩阵（i，j用于表示节点的状态，$P_{ij}=P(x_{t+1}=j|x_t=i)$），${\pi }_t$表示$t$时刻的节点所表示的分布：

13.3.2 Markov Chain性质——平稳分布

Markov Chain性质——平稳分布：节点的分布随时间$t\to \infty$，${\pi }_t$将会收敛于一个固定的分布。

重定义：

• 节点为 { x t } {\lbrace x_t \rbrace} {xt​}， 状态分布为$q^{(t)}(x)$，$Q = \begin{pmatrix} Q_{11} & Q_{12} & \dots & Q_{1K} \ Q_{21} & Q_{22} & \dots & Q_{2K} \ \dots \ Q_{K1} & Q_{K2} & \dots & Q_{KK} \end{pmatrix}_{K \times K} $为转移矩阵（随机矩阵），到$m$时刻收敛（分布不变）：
  

为何要引入平稳分布？

• 从结论出发，平稳分布只与转移矩阵相关（也就是说初始分布可以很简单，简单到可以取样），若平稳分布与目标分布相同的话，我们就可以通过Markov Chain的平稳分布，从简单取样到复杂取样。

收敛证明——证明对于Markov Chain，$\exists m$使得$q^{(m+i)}(x),i\in N$进入平稳分布：

定义行向量：$q^{(t+1)}={\left(q^{(t+1)}(x=1),\dots ,q^{(t+1)}(x=K)\right)}_{1\times K}$

引入转移矩阵的定义可以得到：$q^{(t+1)}(x=j)={\sum }_{i=1}^K{q}^{(t)}(x=i)\cdot Q_{ij}$

可以得到公式：
$$\begin{array}{rl}{q}^{(t+1)}& ={\left(\sum _{i=1}^K{q}^{(t)}(x=i)\cdot Q_{i1},\dots ,\sum _{i=1}^K{q}^{(t)}(x=i)\cdot Q_{iK}\right)}_{1\times K}\\ & =q^{(t)}Q=\cdots =q^{(1)}{Q}^t\end{array}$$
引入随机矩阵的性质（Q满足随机矩阵的性质）：随机矩阵的特征值的绝对值$\le 1$，所以有：
$$Q=A\Lambda A^{-1}=A^{-1}\Lambda A,\Lambda =diag({\lambda }_1,\dots ,{\lambda }_k),|{\lambda }_i|\le 1$$
引入性质继续分解公式可得：
$$\begin{array}{r}{q}^{(t+1)}=q^{(1)}{Q}^t=q^{(1)}{(A\Lambda A^{-1})}^t=q^{(1)}A{\Lambda }^t{A}^{-1}\end{array}$$
假设目前有一个足够大的$m$，我们可以知道：
$$\{\begin{array}{l}{q}^{(m+1)}=q^{(1)}A{\Lambda }^m{A}^{-1}\\ q^{(m+2)}=q^{(1)}A{\Lambda }^{m+1}{A}^{-1}\end{array}$$
若此时$m$足够大，不妨设其中只有${\lambda }_i=1$，可得${\Lambda }^m=diag(0,\dots ,1,\dots ,0)={\Lambda }^{m+1}$​。将该公式代入上文可得
$$q^{(m+1)}=q^{(m+2)}$$
由此得证，当$t>m$时，有$q^{(m+1)}=q^{(m+2)}=q^{(\infty )}$

但这个收敛证明有一个缺点：

• 我们只证明了在足够大时会收敛，但我们不知道具体的收敛时间与条件。

13.3.3 Markov Chain性质——Detailed Balance

Detailed Balance：

• 满足以下等式（表示为双向映射）：
  $$\pi (x)\cdot P(x\mapsto x^{\ast })=\pi (x^{\ast })\cdot P(x^{\ast }\mapsto x)$$

作用：

• Detailed Balance的条件中包含了平稳分布，因为平稳分布的性质表示为： { π k } {\lbrace \pi_k \rbrace} {πk​}是 { x t } {\lbrace x_t \rbrace} {xt​}的平稳分布：
  $$\pi (x^{\ast })=\int \pi (x)\cdot P(x\mapsto x^{\ast })\mathrm{d}x$$

• 要证明Detailed Balance包含平稳分布可以在两边同对$x$​积分：
  $$\begin{array}{rl}\int \pi (x)\cdot P(x\mapsto x^{\ast })\mathrm{d}x& =\int \pi (x^{\ast })\cdot P(x^{\ast }\mapsto x)\mathrm{d}x\\ & =\pi (x^{\ast })\int P(x^{\ast }\mapsto x)\mathrm{d}x\\ & =\pi (x^{\ast })\end{array}$$

13.4 Metropolis Hastings方法

意义：是一种通过接受率满足Detailed Balance条件的方法。

方法：

1. 我们先定义一个随机转移矩阵$Q(Z\mapsto Z^{\ast })$，这个转移矩阵表示条件分布$Q(Z^{\ast }|Z)$，但此时：
   $$p(Z)\cdot Q(Z\mapsto Z^{\ast })\ne p(Z^{\ast })\cdot Q(Z^{\ast }\mapsto Z)$$

2. 定义一个接受率$\alpha (Z,Z^{\ast })=\min (1,\frac{p(Z^{\ast })\cdot Q(Z^{\ast }\mapsto Z)}{p(Z)\cdot Q(Z\mapsto Z^{\ast })})$，使得等式成立：
   $$\begin{array}{rl}p(Z)\cdot Q(Z\mapsto Z^{\ast })\cdot \alpha (Z,Z^{\ast })& =p(Z)\cdot Q(Z\mapsto Z^{\ast })\cdot \min (1,\frac{p(Z^{\ast })\cdot Q(Z^{\ast }\mapsto Z)}{p(Z)\cdot Q(Z\mapsto Z^{\ast })})\\ & =\min (p(Z)\cdot Q(Z\mapsto Z^{\ast }),p(Z^{\ast })\cdot Q(Z^{\ast }\mapsto Z))\\ & =p(Z^{\ast })\cdot Q(Z^{\ast }\mapsto Z)\cdot \min (\frac{p(Z)\cdot Q(Z\mapsto Z^{\ast })}{p(Z^{\ast })\cdot Q(Z^{\ast }\mapsto Z)},1)\\ & =p(Z^{\ast })\cdot Q(Z^{\ast }\mapsto Z)\cdot \alpha (Z^{\ast },Z)\end{array}$$

3. 此时$Q(Z\mapsto Z^{\ast })\cdot \alpha (Z,Z^{\ast })$就相当于是新的转移矩阵。

具体实现：

1. 从均匀分布中取得$u\backsim U(0,1)$
2. 从转移矩阵中进行抽样$Z^{\ast }\backsim Q(Z\mapsto Z^{\ast })$
3. 引入接受率，对样本进行一个筛选：$u\le \alpha ?Z^{(i)}=z^{\ast }:Z^{(i)}=Z^{(i-1)}$

13.5 Gibbs方法

Gibbs实际上是特殊的MH算法，这个时候我们取$\alpha =1$

Gibbs的思想是：

• 固定$p(Z)=p(Z_1,Z_2,\dots ,Z_M)$​中的每个维度，然后进行坐标上升。每个维度就是在条件分布中取样：
  $$Z_i\backsim p(Z_i|Z_{-i})$$

具体实现：

1. 先采样出$Z_1^{(0)},Z_2^{(0)},\dots ,Z_M^{(0)}$

2. 然后不停的上升：
   $$\{\begin{array}{l}{Z}_1^{(t+1)}=p(Z_1|Z_2^{(t)},\dots ,Z_M^{(t)})\\ Z_M^{(t+1)}=p(Z_M|Z_2^{(t+1)},\dots ,Z_{M-1}^{(t+1)})\end{array}$$

3. 这样到最后将维度组合起来就是样本了

为什么Gibbs不用像MH中一样乘$\alpha$，为什么$\alpha =1$：

• 在Gibbs中我们假设$Q(Z\mapsto Z^{\ast })=P(Z_i^{\ast }|Z_i)$，因为我们在求取第$i$项到时候固定住了其他的项，所以有$Z_i=Z_i^{\ast }$，根据这些性质可以得到：
  $$\frac{p(Z)\cdot Q(Z\mapsto Z^{\ast })}{p(Z^{\ast })\cdot Q(Z^{\ast }\mapsto Z)}=\frac{p(Z|Z_i)\cdot p(Z_i)\cdot Q(Z\mapsto Z^{\ast })}{p(Z^{\ast }|Z_i^{\ast })\cdot p(Z_i^{\ast })\cdot Q(Z^{\ast }\mapsto Z)}=1$$

所以我们无需加入$\alpha$

13.6 MCMC总结

整体思路：利用Markov Chain收敛与平稳分布的性质，设计Q使得平稳分布$\approx$目标分布

MCMC有的问题：

• 理论证明了收敛性，但没有给出证明结束时刻m的方法——不能简单判定相邻时刻的分布，如果提早判定了，就有可能无法测出一个较低的高峰：
  

• mixing time过长：$p(x)$过于复杂/高维/高相关性

• 样本间存在一定的相关性。

将MCMC化为一张图就是：