NO.6.tip： 马尔可夫链蒙特卡罗方法

背景

  我们经常需要计算随机变量的某个函数的期望值 $\mathbb{E}[f(\mathbf{X})]$。这需要计算以下积分：

$$\mathbb{E}[f(\mathbf{x})]=\int f(\mathbf{x})p(\mathbf{x})\mathrm{d}\mathbf{x}$$

  其中，$\mathbf{x}\in {\mathbb{R}}^n$，$f:{\mathbb{R}}^n\to {\mathbb{R}}^m$，$p(\mathbf{x})$ 是 $\mathbf{X}$ 的目标分布。在低维度（例如3维）中，我们可以使用数值积分有效地计算上述积分。数值积分（自适应地）计算网格，然后评估网格上每个点的函数。但这无法扩展到更高的维度。

  另一种方法是抽取多个随机样本 ${\mathbf{x}}_n\sim p(\mathbf{x})$，然后计算：

$$\mathbb{E}[f(\mathbf{x})]\approx \frac{1}{N_s}\sum _{n=1}^{N_s}f({\mathbf{x}}_n)$$

  这被称为蒙特卡罗积分。与数值积分相比，蒙特卡罗积分的优势在于，函数只在概率不可忽略的地方进行评估，因此不需要均匀覆盖整个空间。特别是，可以证明，精度原则上与 $\mathbf{x}$ 的维度无关，仅取决于样本数量 $N_s$。关键在于，我们首先需要一种生成样本 ${\mathbf{x}}_n\sim p(\mathbf{x})$ 的方法。

  有几种常见的确定性的采样方式：

逆概率变换：用随机数"反查"概率分布的累积函数，把均匀随机变量变成服从目标分布的样本，从而用样本均值近似任意复杂函数的积分。

拒绝采样：在目标分布难以直接采样时，用一个容易采样的"建议分布"生成候选点，然后按目标分布与建议分布的比值概率接受或拒绝这些点，最终保留下来的样本就服从目标分布，进而用于蒙特卡洛积分。

重要性采样：不直接从目标分布采样，而是从一个更容易采样的"建议分布"中抽取样本，然后给每个样本赋予一个权重（目标密度与建议密度之比）来修正偏差，让积分估计更集中在函数值大的区域，从而降低方差、提高效率。

  确定性方法用规则网格铺满空间，误差 $O(N^{-k/d})$ 随维度指数恶化；蒙特卡洛用随机样本估计期望，误差 $O(N^{-1/2})$ 与维度无关，只取决于样本数N——前者靠"覆盖"换精度，后者靠"概率"换自由。

源与流

  因此，马尔可夫链蒙特卡罗方法就是这样一种流行的从高维分布中采样的方法。马尔可夫链蒙特卡罗方法背后的基本思想是在状态空间 $\mathcal{X}$ 上构建一个马尔可夫链，其平稳分布是感兴趣的目标密度 $p^{\ast }(x)$。（在贝叶斯上下文中，这通常是一个后验，$p^{\ast }(x)\propto p(x\mid \mathcal{D})$，但马尔可夫链蒙特卡罗方法可以应用于从任何类型的分布生成样本。）也就是说，通过从链中抽取（相关的）样本 $x_0,x_1,x_2,\cdots$，我们可以相对于 $p^{\ast }$ 进行蒙特卡罗积分。

  Metropolis-Hastings 算法的基本思想是，在每一步中，我们都建议从当前状态 $x$ 移动到概率为 $q(x^{\prime }\mid x)$ 的新状态 $x^{\prime }$，其中 $q$ 被称为提议分布（也称为核）。用户可以自由使用他们想要的任何类型的提议分布，但须遵守我们在下面解释的一些条件。这使得 Metropolis-Hastings 算法成为一种非常灵活的方法。

  在提议移动到 $x^{\prime }$ 之后，我们根据某种公式决定是接受还是拒绝这一提议，该公式确保在每个状态花费的长期时间与 $p^{\ast }(x)$ 成比例。如果提议被接受，则新状态为 $x^{\prime }$，否则新状态与当前状态 $x$ 相同（即重复样本）。

  如果提议是对称的，那么 $q(x^{\prime }\mid x)=q(x\mid x^{\prime })$，接受概率由以下公式给出：

$$A=\min \left(1,\frac{p^{\ast }(x^{\prime })}{p^{\ast }(x)}\right)$$

  我们看到，如果 $x^{\prime }$ 比 $x$ 更有可能，那么我们肯定会移动到 $x^{\prime }$（因为 $\dfrac{p^{}(x')}{p^{}(x)} > 1$），但如果 $x^{\prime }$ 不太可能，我们仍然可能移动到 $x^{\prime }$，这取决于相对概率。因此，我们不是贪婪地只进入更可能的状态，而是偶尔允许"下坡"进入不太可能的状态。可以证明，该过程可以确保在每个状态 $x$ 下花费的时间等于 $p^{\ast }(x)$。

  如果提议分布是不对称的，那么 $q(x^{\prime }\mid x)\ne q(x\mid x^{\prime })$，我们需要 Hastings 修正（Hastings correction），由以下公式给出：

$$\begin{array}{rl}A& =\min (1,\alpha )\\ \alpha & =\frac{p^{\ast }(x^{\prime })q(x\mid x^{\prime })}{p^{\ast }(x)q(x^{\prime }\mid x)}=\frac{p^{\ast }(x^{\prime })/q(x^{\prime }\mid x)}{p^{\ast }(x)/q(x\mid x^{\prime })}\end{array}$$

  这种修正是为了补偿提议分布本身（而不仅仅是目标分布）可能有利于某些状态的事实。

  Metropolis-Hastings 算法之所以是一种有用算法的一个重要原因是，在评估 $\alpha$ 时，我们只需要知道达到归一化常数的目标密度。特别地，假设 $p^{\ast }(x)=\frac{1}{Z}\tilde{p}(x)$，其中 $\tilde{p}(x)$ 是未归一化分布，$Z$ 是归一化常数。于是：

$$\alpha =\frac{(\tilde{p}(x^{\prime })/Z)q(x\mid x^{\prime })}{(\tilde{p}(x)/Z)q(x^{\prime }\mid x)}$$

  所以 $Z$ 可以消除。因此，即使 $Z$ 未知，我们也可以从 $p^{\ast }$ 中采样。

  Metropolis-Hastings 方法的主要问题是需要选择提议分布，以及接受率可能较低。还有一种经典方法吉布斯采样，该方法利用图模型的条件独立性来自动创建一个好的提议分布，接受概率为1。