一、概述

1.概述

马尔可夫链蒙特卡罗法 (MCMC) 是一种强大的工具，它在处理多元、非标准形式的概率密度函数，以及随机变量各分量不独立的情况时，比拒绝采样法和重要性采样法更具优势。它可以用于从一个概率分布中进行随机抽样，或者计算函数关于该概率分布的数学期望。

假设我们有一个多元随机变量 $x$ ，它满足 $x \in \mathbb{X}$ ，并且其概率密度函数为 $p(x)$ 。我们有一个定义在 $x \in \mathbb{X}$ 的函数 $f(x)$ ，我们的目标是获取概率分布 $p(x)$ 的样本集合，以及计算函数 $f(x)$ 的数学期望 $E_{p(x)}[f(x)]$ 。

MCMC方法可以用来解决这个问题。其基本思想是在随机变量 $x$ 的状态空间 $S$ 上定义一个满足遍历定理的马尔可夫链 $X$ ，使其平稳分布就是我们要抽样的目标分布 $p(x)$ 。然后我们在这个马尔可夫链上进行随机游走，每个时刻都会得到一个样本。根据遍历定理，当时间趋于无穷时，样本的分布会趋于平稳分布，样本的函数均值也会趋近函数的数学期望。

因此，当时间足够长时（即时刻大于某个正整数 $m$ ），在之后的时间（时间小于等于某个正整数 $n ， n>m$ ）里随机游走得到的样本集合 $\left\{x_{m+1}, x_{m+2}, \cdots, x_n\right\}$ 就可以视为目标概率分布的抽样结果，得到的函数均值（遍历均值）就是我们要计算的数学期望值。

\hat{E} f=\frac{1}{n-m} \sum_{i=m+1}^n f\left(x_i\right)

到时刻 $m$ 为止的时间段称为燃烧期，即，从初始状态达平稳分布的时间段。

2. 一些要关注的知识点

由于该马尔可夫链符合遍历定理，随机游走的起点对最终结果没有影响。也就是说，无论起始点在何处，都将会汇集于同一个平稳分布。
MCMC的收敛性判断常常是基于观察的，例如，在马尔可夫链上执行随机游走，并检查遍历平均值是否收敛。具体的手段包括：
- ①间隔一定时间获取一次样本，在获取多个样本后，计算遍历平均值，一旦计算出的平均值稳定，就可以认为马尔可夫链已经收敛。
- ②在马尔可夫链上同时执行多个随机游走，比较各个随机游走的遍历平均值是否相似。
在MCMC中获取的样本序列，相邻的样本点是有关联的，并非独立的。因此，当需要独立样本时，可以在该样本序列中再进行随机抽样，例如，每隔一段时间取一次样本，然后将这样得到的子样本集合当作独立样本集合。
通常情况下，MCMC比拒绝采样更易于实施，因为只需要定义马尔可夫链，而无需定义建议分布。通常情况下，MCMC比拒绝采样更高效，因为没有大量被拒绝的样本，尽管“燃烧期”的成本也是需要舍弃的。

3. 马尔可夫链蒙特卡罗法的核心步骤

①首先，在随机变量 $x$ 的状态空间上设定一个满足遍历定义的马尔可夫链，使其平稳分布为目标分布 $p(x)$ ;
②从状态空间的某一点 $x_0$ 开始，利用设定的马尔可夫链执行随机游走，生成样本序列 $x_0,x_1,\cdots ,x_t,\cdots$ ;
③应用马尔可夫链的遍历定理，确定正整数 $m$ 和 $n$ （ $m<n$ ），得到样本集合 $\left \{x_{m+1},x_{m+2},\cdots ,x_{n}\right \}$ ，然后计算 $f(x)$ 的平均值（遍历平均值）：

$\hat{E}f=\frac{1}{n-m} \sum_{i=m+1}^{n} f(x_{i})$

这里有几个重要问题：

①如何定义马尔可夫链，保证MCMC的条件成立；

②如何确定收敛步数 $m$ ，保证样本抽样的无偏性；

③如何确定迭代步数 $n$ ，保证遍历均值计算的精度。

二、Metropilis-Hastings算法（MH算法）

由于满足Detail Balance的公式即为满足平稳分布（公式迭代为稳态），因此MH算法引入了接受率 $α(z, z^*)$ 的概念。我们提出了一个设想，当未满足Detail Balance等式时，此时状态转移矩阵为Q，正在进行迭代，此时引入接受率，使得等式两边的 Q * 接受率相等，达到最终的稳态。

\begin{aligned} & p(z) \cdot \underbrace{Q\left(z \longmapsto z^*\right) \cdot α\left(z, z^*\right)}_{p\left(z \longrightarrow z^*\right)}=p\left(z^*\right) \cdot \underbrace{Q\left(z^* \rightarrow z\right) \cdot α\left(z^*, z\right)}_{ P\left(z^* \longmapsto z\right)}\\ & p(z) \cdot p\left(z \longrightarrow z^*\right)= p\left(z^*\right) \cdot P\left(z^* \longmapsto z\right)\\ & 定义:α\left(z, z^*\right)=\min \left(1, \frac{p\left(z^*\right) \cdot Q\left(z^* \mapsto z)\right.}{p(z) \cdot Q\left(z \longmapsto z^*)\right.}\right) \\ \end{aligned}

有了接受率，我们可以对当前状态式子进行转化

$\begin{aligned} & p(z) \cdot Q\left(z \mapsto z^*\right) \cdot α\left(z, z^*\right) \quad \\ & =p(z) \cdot Q\left(z \mapsto z^*\right) \cdot \min \left(1, \frac{p\left(z^*\right) \cdot Q\left(z^* \mapsto z^*\right)}{p(z) \cdot Q\left(z \mapsto z^*\right)}\right) \\ & =\min \left(p(z) \cdot Q\left(z \mapsto z^*\right), p\left(z^*\right) \cdot Q\left(z^* \mapsto z\right)\right) \\ & \begin{array}{l} =p\left(z^*\right) Q\left(z^* \mapsto z\right) \cdot \underbrace{\min \left(1, \frac{p(z) \cdot Q\left(z \mapsto z^*\right)}{p\left(z^*\right) \cdot Q\left(z^*\mapsto z\right)}\right.}) \\ =p\left(z^*\right) \cdot Q\left(z^* \mapsto z\right) \cdot α\left(z^*, z\right) \end{array} \\ & \end{aligned}$
验证构造的接受率能等于原式，形成闭解

Metropolis-Houstings:

$u \sim U(0,1)$ ，U为均匀分布的意思

$z^* \sim Q\left(z \mid z^{(i-1)}\right)$ (Q可随便取随机矩阵)

$α=\min \left(1 . \frac{p\left(z^*\right) \cdot Q\left(z^* \mapsto z\right)}{p(z) \cdot Q\left(z \rightarrow z^*\right)}\right)$ if $u \leqslant α, z^{(i)}=z^*$ else $z^{(i)}=z^{(i-1)}$

上述if-else的意思是，若随机采取的样本u小于α时，就将该 $z^*$ 就被接受，作为第i个样本，否则不接受该样本，使用上一样本重新刷新u到下一状态

假如随机执行迭代N次（生成n个u）可得到 $z^{(1)}, z^{(2)}, \cdots, z^{(N)}$ 个样本点，这样我们可以认为这些样本点来自于后验分布 $P(Z)$ ，这样的话可以把n个样本点带入，求得我们最终的目标

\begin{aligned} p(z) \rightarrow & E_{p(z)}[f(z)]=\int_z f(z) p(z) d z \\ & \approx \frac{1}{N} \sum_{i=1}^N f(z^{(i)}) \end{aligned}

关于 $p(z)$ 的补充：

$p(z)=\frac{\hat{p}(z)}{z_p}$ ，其中 $z_p$ 为归一化因子，不需要求（α式子中上下抵消），而 $\hat{p}(z)$ 为估计，估计通常都能求，一般就是likelihood * prior即可，也就是我们并不是直接使用 $p_{(z)}，p_{(z^*)}$ ，所以说，真正的式子应为 $\frac{\hat{p}\left(z^*\right) \cdot Q\left(z^* \mapsto z\right)}{\hat{p}(z) \cdot Q\left(z \rightarrow z^*\right)}$

即，样本点 $z^{(i)}$ ~ $\hat{p}(z)$ ~ $p_{(z)}$

三、吉布斯抽样

吉布斯抽样可以认为是MH算法的特殊情况，但是更容易实现，因此被广泛使用。

与坐标上升法类似，求一个参数时确定别的参数为定值。

Gibbs $\quad p(z)=p\left(z, z_2, \cdots, z_m\right)$

$z_i \sim p\left(z_i \mid z_{-i}\right) \\ \hookrightarrow z_1, z_2, \cdots, z_{i-1}, z_{i+1}, \cdots z_M \\ 设样本一共有三个\\$

$\begin{aligned} & p(z)=p\left(z_1, z_3, z_3\right) \\ & i: \quad z_1^{(0)}, z_2^{(0)}, z_3^{(0)} \\ & \text { ii: } t+1 \text {， } \\ & z_1^{(t+1)} \sim p\left(z_1 \mid z_2^{(t)} z_3^{(t)}\right) \\ & z_2^{(t+1)} \sim p\left(z_2 \mid z_1^{(t+1)}, z_3^{(t)}\right) \\ & z_3^{(t+1)} \sim p\left(z_3 \mid z_1^{(t+1)} , z_2^{(t+1)}\right) \\ & 注意观察上述三个式子的规律，先得到z_1^{(t+1)}，就将其运用于z_2^{(t+1)}...\\ & \frac{p\left(z^*\right) \cdot Q\left(z^* \mapsto z\right)}{p(z) \cdot Q\left(z \mapsto z^*\right)} \\ & =\frac{p\left(z_i^* \mid z_{-i}^*\right) \cdot p\left(z_{-i}^*\right) p\left(z_i \mid z_{-i}^*\right)}{p\left(z_i \mid z_{-i}\right) \cdot p\left(z_{-i}\right) \cdot p\left(z_i^* \mid z_{-i}\right)}=1 \\ & \end{aligned}$
采样的时候其余纬度固定，所以 $z^*_{-i} = z_{-i}$ ，这也是为什么上式能等于1的原因

显然它的效率比传统的MH要高，因为接受率永远是1，吉布斯抽样对每次抽样的结果都接受，没有拒绝，这一点和一般的MH算法不同。

对比吉布斯抽样和MH：

单分量MH算法与吉布斯抽样的主要差别在于，前者的抽样可能在样本点之间跳动，但有可能在某些样本点上停留（因为采样被拒绝）；而后者的抽样点则会在样本点之间连续移动。

吉布斯抽样适合于条件概率分布可以容易抽样的情况，而单分量MH算法则适合于条件概率分布抽样困难的情况，此时使用容易抽样的条件分布作为建议分布。

四、采样思路以及面临困难

思路：利用 Markov Chain 收敛于平稳分布。设计Q，使得平稳分布(q(x)) ≈ 目标分布(p(x))

设离散的马尔科夫链：

状态空间: $\{ 1,2,3,...,k \}$

状态转移矩阵： $Q = [Q_{ij}]_{k×k}$ (连续链，即为转移概率密度)

给定状态转移矩阵

Q=\begin{pmatrix} Q_{11}&Q_{12}&\cdots&Q_{1K}\\ \vdots&\vdots&\vdots&\vdots\\ Q_{i1}&Q_{i2}&\cdots&Q_{iK}\\ \vdots&\vdots&\vdots&\vdots\\ Q_{k1}&Q_{k2}&\cdots&Q_{KK}\end{pmatrix}_{k×k}

那么 $p(x^{(2)}|x^{(1)}=i)$ 即为 $Q_{i1}\;\;Q_{i2}\cdots Q_{iK}$

问题：

① 理论只保证收敛性，但无法知道何时收敛
② mixing time 过长 --》 p(x)太复杂 --》高维特征之间的相关性导致（例如GMM，可能跨不过山谷，而取样集中于某一高斯分布）
③ 样本本身是取自各个m时刻后的状态，本身可能就一定相关性，毕竟每个状态之间都有转移矩阵。

“开启掘金成长之旅！这是我参与「掘金日新计划 · 2 月更文挑战」的第 12 天，点击查看活动详情”

15.马尔可夫链蒙特卡罗法（MCMC）：一种强大的统计工具

一、概述

二、Metropilis-Hastings算法（MH算法）

三、吉布斯抽样

四、采样思路以及面临困难