马尔科夫蒙特卡洛（MCMC）模拟

问题：找后验分布

假设我们有一组随机变量 $\bm{\theta}$，以及观测数据 $\bm{y}$，我们想要找到参数的后验分布 $p(\bm{\theta \ | \ y})$。根据贝叶斯定理，可以将后验分布表示为：

$$p (\bm{\theta \ | \ y}) = \frac{p (\bm{y \ | \ \theta}) \times p(\bm{\theta})}{\int p (\bm{y \ | \ \theta}) \times p(\bm{\theta}) \ d\bm{\theta}}$$

其中，$p (\bm{y \ | \ \theta})$ 是假设$\bm{\theta}$下的观测到$\bm{y}$的概率，可以根据模型来计算；$p(\bm{\theta})$ 是先验分布，一般是假设的。

在实际应用中，分母的计算通常非常困难，因为它涉及到对参数空间进行积分，所以分母很难算，但是分母是一个常数，即使不管它也不影响后验分布的形状，所以可以简化表示为：

$$p (\bm{\theta \ | \ y}) \propto p (\bm{y \ | \ \theta}) \times p(\bm{\theta})$$

分子是可以计算的，因为先验分布$p(\bm{\theta})$的概率密度函数是已知的，而$p (\bm{y \ | \ \theta})$也可以通过模型来计算。

所以理论上只需要对参数空间进行均匀采样，比如10均分，得到一系列参数$\bm{\theta}$的样本，然后根据这些样本计算分子，就能得到后验分布的近似值，也能找到后验分布的峰值。但是在高维空间中，均匀采样的效率非常低，比如如果有10个参数，每个参数均匀采样10个点，那么就需要$10^{10}$次采样才能覆盖整个参数空间，这在计算上是不可行的。

马尔科夫链蒙特卡洛（MCMC）模拟

除了均匀采样，还可以使用马尔科夫链蒙特卡洛（MCMC）模拟直接从后验分布中采样，只需要知道后验分布的比例关系（即分子）。

马尔科夫链是指一组随机状态，其中每个状态只依赖于前一个状态，所以一个马尔科夫链的状态一直转移，样本多了之后，就会收敛到一个稳定的分布 $\pi(\bm{\theta})$，$\pi(\bm{\theta})$表示$\bm{\theta}$的概率密度函数。

MCMC模拟就是通过构造一个马尔科夫链来生成一系列参数$\bm{\theta}$的样本，当样本数量足够多时，这些样本就会近似于后验分布 $p (\bm{\theta \ | \ y})$。关键是如何构造这个马尔科夫链，使得它的稳定分布就是我们想要的后验分布。

Metropolis-Hastings 算法

Metropolis-Hastings 算法是这样实现MCMC模拟的：

1. 从一个当前状态开始 $\bm{\theta}_{\mathrm{current}}$。
2. 提议一个新的状态 $\bm{\theta}_{\mathrm{new}}$，这里提议分布可以是任意的，表示为 $q(\bm{\theta}_{\mathrm{new}} | \bm{\theta}_{\mathrm{current}})$。
3. 计算接受概率 $\alpha$：

   $$\alpha = \frac{p (\bm{y \ | \ \theta}_{\mathrm{new}}) \times p(\bm{\theta}_{\mathrm{new}}) \times q(\bm{\theta}_{\mathrm{current}} \ | \ \bm{\theta}_{\mathrm{new}})}{p (\bm{y \ | \ \theta}_{\mathrm{current}}) \times p(\bm{\theta}_{\mathrm{current}}) \times q(\bm{\theta}_{\mathrm{new}} \ | \ \bm{\theta}_{\mathrm{current}})}$$

   这里的接受概率是一个精心设计的比例关系，确保了马尔科夫链的稳定分布是我们想要的后验分布，后面证明这个接受概率的设计是合理的。
4. 如果 $\alpha \geq 1$，则接受并转移到新的状态 $\bm{\theta}_{\mathrm{new}}$。
5. 如果 $\alpha < 1$，则以概率 $\alpha$ 接受新的状态 $\bm{\theta}_{\mathrm{new}}$，否则保持当前状态 $\bm{\theta}_{\mathrm{current}}$，但是无论如何都将当前状态添加到样本中。
6. 重复直到收集到足够的样本。

结果生成了一个马尔科夫链的样本，也就是一个参数$\bm{\theta}$的样本集合，这些样本会近似于后验分布 $p (\bm{\theta \ | \ y})$。

证明为什么上面的马尔科夫链会收敛到后验分布

定性分析

接受概率 $\alpha$ 是一个比值。

分子是 $p (\bm{y \ | \ \theta}_{\mathrm{new}}) \times p(\bm{\theta}_{\mathrm{new}}) \times q(\bm{\theta}_{\mathrm{current}} \ | \ \bm{\theta}_{\mathrm{new}})$，前面两项其实就是是后验分布 $\pi(\bm{\theta}_{\mathrm{new}}) = p (\bm{y \ | \ \theta}_{\mathrm{new}}) \times p(\bm{\theta}_{\mathrm{new}})$。第三项当提议分布是均匀的时候为常数。所以分子表示了新状态 $\bm{\theta}_{\mathrm{new}}$ 的后验概率密度，分母表示了当前状态 $\bm{\theta}_{\mathrm{current}}$ 的后验概率密度。

如果新状态的后验概率密度更高，那么接受概率 $\alpha$ 就大于1，必然接受新状态；如果新状态的后验概率密度更低，那么接受概率 $\alpha$ 就小于1，以一定概率接受新状态，这样就保证了马尔科夫链会倾向于向后验分布的高概率区域移动。所以最终高概率区域的样本会更多，低概率区域的样本会更少，这样就能近似于后验分布。

证明

严格一点的证明是需要利用细致平衡（detailed balance）条件来证明的。

对于一个马尔科夫链，转移概率为 $T(\bm{\theta}_a \to \bm{\theta}_b)$，如果任意两个状态 $\bm{\theta}_a$ 和 $\bm{\theta}_b$，满足：

$$\pi(\bm{\theta}_a) T(\bm{\theta}_a \to \bm{\theta}_b) = \pi(\bm{\theta}_b) T(\bm{\theta}_b \to \bm{\theta}_a)$$

那么分布 $\pi(\bm{\theta})$ 就是这个马尔科夫链的稳定分布。

上面的细致平衡条件可以理解为，左边是状态 $\bm{\theta}_a$ 附近单位长度的粒子流向状态 $\bm{\theta}_b$ 的数量，右边是状态 $\bm{\theta}_b$ 附近单位长度的粒子流向状态 $\bm{\theta}_a$ 的数量，如果这两个概率相等，那么就说明在这个马尔科夫链中，状态 $\bm{\theta}_a$ 和 $\bm{\theta}_b$ 之间的粒子流是平衡的。

后面只需要证明 Metropolis-Hastings 算法所构造的转移概率满足细致平衡条件，就能证明这个马尔科夫链的稳定分布是我们想要的后验分布。

对于任意两个状态 $\bm{\theta}_a$ 和 $\bm{\theta}_b$，Metropolis-Hastings 算法的转移概率有两种情况：

$$T(\bm{\theta}_a \to \bm{\theta}_b) = \begin{cases} q(\bm{\theta}_b \mid \bm{\theta}_a) \cdot \alpha, & \alpha < 1 \\ q(\bm{\theta}_b \mid \bm{\theta}_a), & \alpha \geq 1 \end{cases}$$

其中

$$\alpha = \frac{p (\bm{y \mid \theta}_b) \, p(\bm{\theta}_b) \, q(\bm{\theta}_a \mid \bm{\theta}_b)}{p (\bm{y \mid \theta}_a) \, p(\bm{\theta}_a) \, q(\bm{\theta}_b \mid \bm{\theta}_a)}$$

对于 $\alpha < 1$ 的情况，$T(\bm{\theta}_a \to \bm{\theta}_b)$取上面的表达式，但是$T(\bm{\theta}_b \to \bm{\theta}_a)$刚好相反为$q(\bm{\theta}_a \mid \bm{\theta}_b)$，代入细致平衡条件：

$$\pi(\bm{\theta}_a) q(\bm{\theta}_b \mid \bm{\theta}_a) \cdot \alpha = \pi(\bm{\theta}_b) q(\bm{\theta}_a \mid \bm{\theta}_b)$$

里面的 $\pi(\bm{\theta}_a) = p (\bm{y \mid \theta}_a) \, p(\bm{\theta}_a)$，$\pi(\bm{\theta}_b) = p (\bm{y \mid \theta}_b) \, p(\bm{\theta}_b)$，代入发现 $\alpha$ 的定义正好满足细致平衡条件。

对于 $\alpha \geq 1$ 的情况，也是类似的。