第2章 · 蒙特卡洛基础#

蒙特卡洛方法建立在概率论与数理统计（probability and statistics）的基础上。本章首先介绍理解从概率分布中采样这一概念所必需的概率论与数理统计基本概念，然后讨论几种实用的技术，用于从仅依赖于少量随机变量的分布函数中进行采样。依赖于大量随机变量的函数需要不同的技术，对这类函数我们将讨论马尔可夫链采样（Markov chain sampling）的使用。

2.1 概率论基本概念#

作为一种随机过程，蒙特卡洛模拟生成一组随机事件。有时事件生成的顺序很重要，有时则不。蒙特卡洛采样理论——即我们获取随机事件的过程——基于概率论。我们先讨论概率论中的几个基本概念。

在概率论中，一个真实或虚拟实验的所有可能结果 $\{\chi_{1} , \chi_{2} , . . . , \chi_{n} , . . . \}$ 构成了一个样本空间（sample space）。例如，如果我们的实验是掷一对骰子，样本空间中就有36种结果。事件（event）是该空间中满足某种条件的一个或多个结果组成的集合。事件A的概率（probability）是分配给该事件的一个数字，需满足三条公理：（1）$0 \leq P ( A ) \leq 1$；（2）如果事件包含样本空间中的所有可能结果，则P(A) = 1；（3）如果事件A分解为两个没有共同结果的事件B和C，则 $P ( A ) = P ( B ) + P ( C )$。在我们的掷骰子实验中，如果事件A是那些总和等于3的结果，它包含两个结果，因此 $P ( A ) = {\frac{1} {18}}$。如果事件A是那些总和为奇数的结果，则 $\begin{array} {r} {P ( A ) = \frac{1} {2}} \end{array}$。

这些公理的一个推论是：概率将事件映射到区间[0, 1]，并且对于一组覆盖样本空间的互斥事件 $\{A_{1} , A_{2} , \ldots \}$，其概率之和等于1。两个事件 $A_{i}$ 和 $A_{j}$ 互斥（mutually exclusive），即没有共同结果，当且仅当 $A_{i}$ 的发生意味着 $A_{j}$ 不发生，反之亦然。由于总和为1，至少有一个事件是可能的。

随机变量（random variable）X将结果映射为一个实数。事件

分别是所有结果 $\chi$ 的集合，使得 $X ( \chi )$ 小于等于 $x$、介于 $x_{1}$ 和 $x_{2}$ 之间、等于 $x_{0}$。通常我们在随机变量层面上工作，即我们的实验生成带有数值的事件。在下文中，我们将概率 $P ( X ( \chi ) = x )$ 写作 $P ( x )$ 或 $P ( X )$，即随机变量及其值可互换使用。此外，我们常常不区分结果（基本事件）与事件。

随机变量X的累积分布函数（cumulative distribution function）$F_{X} ( x )$ 定义为

更一般地，

$F_{X} ( x )$ 是定义在整个实轴上的x的正、单调非递减函数，具有性质 $F_{X} ( - \infty ) = 0$ 和 $F_{X} ( \infty ) = 1$。如果 $F_{X} ( x )$ 处处可微，则随机变量是连续的。此时存在一个概率密度函数（probability density）$f_{X} ( x )$，满足

该表达式说明，$f_{X} ( x ) d x$ 是 $x$ 落在区间 $[ x , x + d x ]$ 内的概率。按照惯例，$F_{X} ( x )$ 是右连续的，即 $F_{X} ( x ) =$ $\mathrm{lim}_{\epsilon \to 0} F_{X} ( x + \epsilon )$。如果 $F_{X} ( x )$ 是阶梯连续的，在 $x_{1} , x_{2} , \dotsc$ 处有跳跃 $f_{1} , f_{2} , \ldots$，则随机变量是离散的。对于离散情形，我们常将 $f_{X} ( x_{i} )$ 简写为 $f_{i}$。我们不显式考虑混合型随机变量的分布，但处理起来相对直接。

有时候在不同随机变量之间进行变换是很有用的。例如，考虑 $Y = y ( X )$，其中 $y ( x )$ 是 $x$ 的非递减函数。如果我们已知 $F_{X} ( x )$ 和 $f_{X} ( x )$，那么 $F_{Y} ( y )$ 和 $f_{Y} ( y ) {\smash [ t ] {\mathstrut}}$ 是多少？由于当 $X \leq x$ 时 $y ( X ) \leq y ( x )$，反之亦然，因此有

即 $F_{Y} ( y ) = F_{X} ( x )$。通过微分，我们得到

因为 $y ( x )$ 是非递减的，所以 $y$ 对 $x$ 的导数为正，方程两边均为正。

如果变量变换是非递增的，很容易证明

且

更一般地，我们写作

或者对于逆映射 $\begin{array} {r} {f_{Y} ( y ) = f_{X} ( x ) \left| \frac{d x} {d v} \right|} \end{array}$

如果我们有两个随机变量 X 和 Y，我们将两者同时发生的概率表示为 $P ( X , Y )$。显然，$P ( Y , X ) = P ( X , Y )$。这个二元函数称为 X 和 Y 的联合概率（joint probability），它与 P(X) 和 $P ( Y )$ 的关系为

如果 $P ( X , Y ) = P ( X ) P ( Y )$，则称这些随机变量是统计独立的（statistically independent）。

在已知事件 Y 发生的条件下，X 的条件概率（conditional probability）$P ( X | Y )$ 与联合概率的关系为

类似地，已知事件 X 发生条件下 Y 的条件概率为

通过比较这两个方程，我们得到

因而有

这被称为贝叶斯定理（Bayes's theorem）。注意

统计独立性意味着 $P ( X | Y ) = P ( X )$，前提是 $P ( Y ) > 0$。

二元累积分布（bivariate cumulative distribution）定义为

与 (2.2) 对应的结果为

相关的密度函数为

对于多于两个随机变量的分布和密度函数，可以类似地定义。

如果我们将两个随机变量变换为新的随机变量 U 和 V，

且逆映射为

那么它们的密度函数之间的关系为

注意雅可比行列式的绝对值。

在多变量分布中，经常需要对一个或多个随机变量进行积分，其结果是一个新的分布，称为边际概率密度（marginal probability density）。例如，设 $f_{X Y} ( x , y )$ 是两个随机变量 X 和 Y 的联合分布，我们关心 $a < X < b$ 的概率。此事件仅在 $a < X < b$ 且 $- \infty < Y < \infty$ 时发生。对于连续和离散情况：

积分或求和的结果是仅关于 $x$ 的函数：

这个函数 $f_{X} ( x )$ 就是 $X$ 的边际概率分布（marginal probability distribution）。

利用边际分布（marginal distribution）的概念和条件概率分布函数（conditional probability distribution function）的多元定义：

我们可以推导出一些有趣的关系。其中之一是查普曼-科尔莫戈罗夫方程（Chapman-Kolmogoroff equation）：

为了推导它，我们使用条件密度（conditional density）的定义，将上式左边写为：

再次使用这个定义并调用边际分布（marginal distribution）的定义，我们得到：

我们可以用类似的推理方法证明诸如：

我们还注意到，重复应用多元条件概率（multivariate conditional probability）的定义可以得到条件概率链式法则（conditional probability chain rule）：

2.2 随机抽样#

我们需要的最基本的蒙特卡洛概念是随机抽样（random sampling）。随机抽样，或简称为抽样，连接了概率论和统计学。概率概念中隐含着存在某种实验，可以从一组可能结果中产生一个随机变量的实现。

如果执行这个实验 $N$ 次，在 $N_{x}$ 次中产生随机变量 $X$ 的值 $x$，那么比值 $N_{x} / N$ 就是 $P ( x )$ 的一个经验度量。对一个概率分布进行抽样，意味着我们以与其成正比的频率生成事件。

随机抽样的另一种方法是同时考虑 $N$ 个实验（一个系综）。与第 $i$ 个实验相关联的是随机变量 $X_{i}$，它与随机变量 $X$ 相同，并且根据假设与 $X$ 具有相同的分布；即 $P_{i} ( X_{i} ) = P ( X )$ 。这 $N$ 个随机变量的联合分布是 $P ( X_{1} , X_{2} , \ldots , X_{N} )$ 。随机变量的统计独立性（statistical independence）意味着：

由N次实验构成的系综因此产生一组结果 $\chi = \{\chi_{1} , \chi_{2} , . . . ,$ ${\chi}_{N} \big \}$，其中第i个随机变量取值 $X_{i} ( \chi ) ~ = ~ X ( \chi_{i} )$ 。独立性与同分布的结合是随机样本（random sample）的典型特征。

对概率分布进行抽样有多种方法，它们几乎都假设存在一个随机数生成器，该生成器是一种数值程序，能产生在区间[0, 1]上均匀分布的随机变量。¹ 换句话说，它在此区间上对 $u ( x ) d x \ : = \ : d x$ 进行抽样。² 实际上，生成的数字是伪随机数（quasi-random）而非真正的随机数：如果长时间对生成器进行采样，这些数字会重复出现。对大多数应用而言，有效序列的长度是足够的。

在图2.1中，我们示意性地展示了典型随机数生成器产生的直方图类型。该直方图在统计涨落范围内是平坦的。如果随机变量在[0, 1]上均匀分布，那么它们的数对应在单位正方形上均匀分布；三元组在单位立方体上均匀分布，依此类推。一般而言，随机数生成器的问题不在于无法产生足够平坦的直方图，而在于能否产生独立的随机数。在数对、三元组或更高阶元组之间可能存在相关性。关于什么是好的随机数生成器以及如何构建它，已有大量文献论述。我们不讨论这些问题，只是假设我们能够获得一个优秀的随机数生成器。如今，现代编译器附带的供应商提供的生成器以及近期教材中介绍的其他生成器已经足够使用。对于大规模严谨模拟，可能需要其他选择。

图2.1 随机数生成器返回的均匀分布示意直方图。生成器被采样50,000次，结果累积到30个箱中。每个箱的计数除以50,000/30（即每个箱的预期计数）。离散概率的片段叠加在均匀分布上。

直接抽样法（Direct sampling）和马尔可夫链方法（Markov chain method）是与量子蒙特卡罗算法相关的两类抽样方法。现在我们讨论每类方法中用于离散和连续分布函数抽样的各种技术。

2.3 直接抽样方法#

2.3.1 离散分布#

直接抽样方法中最简单的例子是从N个事件的离散概率函数 $f_{i}$ 进行抽样的标准程序。为了阐述这一方法，我们假设 $N = 3$，且 $f_{1} = 0.3$，$f_{2} = 0.5$，$f_{3} = 0.2$。接着，我们将这些概率布置在如图2.1所示的均匀分布上。它们将[0, 1]区间划分为三个区域。从该图可以看出，$f_{1}$ 覆盖了随机数生成器产生的大约30%的事件，$f_{2}$ 覆盖50%，$f_{3}$ 覆盖剩余的20%。因此，如果我们从生成器中抽取一个均匀分布的随机数 $\zeta$，我们观察到：如果 $\zeta$ 小于或等于 $f_{1}$，它正确地对事件1进行抽样；如果它大于 $f_{1}$ 但小于或等于 $f_{1} + f_{2}$，则对事件2抽样；如果它大于 $f_{1} + f_{2}$ 但小于或等于 $f_{1} + f_{2} + f_{3} = 1$，则对事件3抽样。这一观察转化为一个简单的策略：利用累积概率函数（cumulative probability function）$F_{i}$ 从离散分布中抽样。我们在算法1中详细描述这一策略。

算法1 通过累积分布从离散分布函数中抽样。 输入：概率向量 $f$。 $F ( 0 ) \gets 0$； 对于 i = 1 到 N 执行 $F ( i ) f ( i ) + F ( i - 1 )$； 创建累积分布函数 结束循环 生成一个均匀随机数 $\zeta \in [ 0 , 1 ]$ $k \gets 0$； 重复 $k \gets k + 1$ 直到 $\zeta \leq F ( k )$ 返回 $k .$

更正式地说，对于均匀随机变量，当 $0 \leq x \leq 1$ 时有 $u ( x ) = f_{X} ( x ) = 1$，因此其累积分布函数（cumulative distribution function）为 $F_{X} ( x ) = x .$。我们还有

这表明$\zeta$落在[0, 1]区间$[x_{1}, x_{2}]$内的概率与该区间的长度成正比。现在，如果我们有一组离散事件，其概率为$f_{i} , i = 1 , \ldots , n$，并且我们希望随机抽取其中一个，我们可以将区间[0, 1]划分为长度为$f_{i}$的子区间。$\zeta$落在哪个子区间，就选择对应的事件。我们通过找到满足下式的k来选择事件k：

如果离散概率只有少数几个元素，那么如算法1所述的简单线性搜索算法就足够了。对于元素数量较多的情况，使用二分查找可以提高效率。如果概率向量在每次抽样时都会变化，可以轻松修改算法1，以省去创建和存储累积分布的需要；参见算法2。如果概率向量非常大且变化不频繁（甚至不变），则可能首选其他抽样算法，如分割法和别名法（参见附录A中的Walker别名法）。

算法1提供了一种从玻尔兹曼概率(Boltzmann probability) $p ( C )$ (1.8)中抽样伊辛构型(Ising configuration) $C_{i}$的方法。在构建好$p ( C_{i} )$之后，我们可以使用(2.10)来抽样$C_{i}$。正如我们之前提到的（第1.3节），当晶格尺寸变大时，构建所有$p ( C_{i} )$的变得不切实际。对于大型晶格，从玻尔兹曼密度中抽样构型通常可以通过马尔可夫链(Markov chain)高效地完成，我们将在下一节讨论这种方法。

算法2 抽样一个离散分布函数。 输入：概率向量 f。 生成一个均匀随机数 $\zeta \in [ 0 , 1 ]$ $k \gets 0$； 重复 $k \gets k + 1$； $\zeta \gets \zeta - f ( k )$； 直到 $\zeta < 0$ 返回 $k .$

一些物理问题具有非常特定的离散分布，例如泊松分布(Poisson distribution)和二项分布(binomial distribution)，这些分布需要进行抽样。累积概率方法(2.10)可用于从中抽样。对于这些以及许多其他离散分布，利用其特性通常更高效（Everett and Cashwell, 1983; Fishman, 1996）。例如，虽然我们可以通过下式抽样泊松概率（图2.2）：

通过

图2.2 泊松分布(2.11)的(a)离散概率和(b)累积概率。这里$t = 1$

一种更高效的方法是

其中$\zeta_{i}$是来自随机数生成器的连续抽取结果。我们将后一种方法的推导留作练习。同时，很容易证明泊松分布(2.11)的均值等于t，方差也等于t。

方程(2.13)定义了一种确定整数k的方法。然而，这个整数是t的函数。在泊松分布的许多应用中，t定义了某个区间的长度，例如(0, t)，而我们感兴趣的不仅是该区间内发生的事件数量$k ( t )$，还包括与这些事件相关的一个递增的点序列$0 < t_{1} < t_{2} \ldots < t_{k} < t$。生成这些$k ( t )$个点的随机过程称为泊松过程(Poisson process)。我们将在第5.2.3节讨论这个过程。

2.3.2 连续分布#

当我们从离散分布的抽样转向连续分布时，我们的注意力也从概率函数转向概率密度(probability density) $f ( x )$。我们可以类比离散分布的程序来抽样$f (x)$。在抽取随机数后，我们求解

对于 $x$。该表达式的内容是 $f ( x ) d x = d \zeta$；也就是说，$x$ 在区间 $[ x , x + d x ]$ 上的概率对应于随机数在区间 $[ \zeta , \zeta + d \zeta ]$ 上的概率。这是连续分布的基本抽样概念。

凭直觉可能足以证明此过程以及离散分布过程的正确性。从数学上讲，两者都需要从以下观察出发：随机变量的函数是另一个随机变量，且累积分布的范围是 [0, 1]，这与我们的均匀随机变量 $\zeta$ 的定义域相同。接下来，我们考虑非递减函数 $F ( x )$（例如累积分布）的广义逆的性质。

非递减函数的广义逆（图 2.3）定义为

图 2.3 广义逆定义的图形表示。

一般来说，该逆函数满足

因此，满足 $F^{-} ( y ) \leq x$ 的点集 $\{( x , y ) \}$ 与满足 $F ( x ) \geq y$ 的点集相同。将 $y$ 与 $\zeta$ 等同，并将 $F ( x )$ 与随机变量 $X$ 的累积分布等同，我们注意到如果 $y$ 在 [0, 1] 中，那么 $x$ 在 $[ F^{-} ( 0 ) , F^{-} ( 1 ) ]$ 中。然后利用累积分布的定义 (2.1)，我们可以表述

因此，要生成一个具有分布 $F ( X )$ 的随机变量 $X$，只需生成一个均匀随机变量 $\zeta$，然后进行变换 $x = F^{-} ( \zeta )$ 即可。

从目前来看，(2.14) 很难直接解出 $x$。如果 $F ( x )$ 是 $f ( x )$ 积分的显式函数表示，那么问题（正如刚才论证的）就简化为求解方程 $\zeta = F ( x )$ 中的 $x$。这也可能很困难，但如果 $F ( x )$ 的解析逆已知，那么

(2.14) 最简单的应用是从区间 $[ a , b ]$ 上的均匀密度进行抽样：

解出 $x$，我们得到

这个结果我们很容易就能猜到。

(2.15) 的一个更有用的应用是从定义在区间 $\lbrack 0 , \infty )$ 上的指数密度进行抽样。我们对其积分得到：

然后解出 $x$：

如果 $\zeta$ 是 (0, 1) 上的随机变量，那么 $1 - \zeta$ 也是。因此，结果是 $x = - \ln \zeta / \lambda$。

我们有时可以通过变量替换来对其他密度进行抽样。让我们使用这种技巧和一个巧妙方法推导出用于从高斯分布（Gaussian distribution）抽样的 Box-Muller 方法：

这个巧妙方法在于观察到从两个高斯分布的乘积中抽样更容易：

通过变换到极坐标：

我们可以写出：

现在我们需要抽样两个概率密度，因此需要两个随机数。对于依赖于 $r$ 的分布，我们有：

由此可得 $r = [ - 2 \log \zeta_{1} ]^{1 / 2}$。对于 $\theta$ 分布，我们需要在 [0, 2π] 上抽样一个均匀分布，得到：

根据两个随机数，得到两个独立的高斯随机变量：

在附录 B 中，我们将讨论拒绝法（rejection method），这是一种从少数几个随机变量的连续分布中采样的通用技术。该技术的一个重要特点是，被采样的分布无需是归一化的。在下一节中，我们将开始讨论马尔可夫链（Markov chain）以及生成它们的常用方法。使用这些链几乎是我们对非常高维的离散和连续分布进行采样的通用方式。

2.4 马尔可夫链蒙特卡洛（Markov chain Monte Carlo）#

在定义马尔可夫链之前，我们先做一些符号上的调整和几点说明。在我们关于马尔可夫链以及后续更一般的蒙特卡洛方法的讨论中，将随机变量 X 最好视为一个“随机向量”$X = ( Y_{1} , Y_{2} , \ldots , Y_{n} , . ~ . ~ . )$，其分量 $Y_{i}$ 是随机变量。这种理解通过消除冗长的下标串，简化了定义多元分布采样过程的符号。这些 $Y_{i}$ 是独立的随机变量，因为我们可以独立地改变每个变量的值，而不受其他变量的影响。随着 $Y_{i}$ 的变化，它们会导致 X 在问题的样本空间内变化。然而，马尔可夫链生成的 $Y_{i}$ 序列之间不一定具有统计独立性。

让我们回想一下第 1 章中对伊辛模型（Ising model）的讨论。在那里，我们将每个伊辛自旋视为一个随机变量，取值为正一或负一。一个构型或状态对应于 $2^{N}$ 个可能结果样本空间中的一个结果 $( s_{1} , s_{2} , \ldots , s_{N} )$。伊辛模型的能量 $E ( C )$ 是一个将构型映射到实数的随机变量。在这个例子中，随机变量 X 代表一个构型 $C$，而 $Y_{i}$ 代表单个自旋，例如 $s_{i}$。用伊辛变量的各个值来表示能量（以及其他可观测量）是很方便的，因此将 C 视为集合 $( s_{1} , s_{2} , \ldots , s_{N} )$。此外，在马尔可夫链中，从一个构型到下一个构型的转移最容易通过一次改变构型的一个分量来实现。

马尔可夫链的定义和性质适用于经典和量子蒙特卡洛算法。在量子蒙特卡洛模拟中，我们可能需要区分系统的状态和蒙特卡洛构型。在量子力学中，我们处理哈密顿算符 H 和状态矢量 $| \psi \rangle$。为了数值地处理它们，我们需要参考某个基矢。给定一组“构型态” $\{| C \rangle \}$，我们可以计算矩阵元 $\langle C^{\prime} | H | C \rangle$ 和波函数 $\psi ( C ) = \langle C | \psi \rangle$，这些是经典对象，因此适合于我们接下来要讨论的采样类型。

在前几节中，我们区分了离散随机变量（具有有限个值）和连续随机变量（具有无穷多个值）。继续这样区分会变得繁琐。在接下来的几节中，符号对于离散随机变量最合适，但在大多数情况下，它可以立即推广到连续情况。然而，我们引用或证明的一些结果，只有在有限、离散时间的马尔可夫链中才是明显正确的。

2.4.1 马尔可夫链（Markov chain）#

马尔可夫链是一种生成随机变量值序列 $x_{1} , x_{2} , \ldots , x_{j} , . . .$ 的过程。如果相关随机变量 $X_{j}$ 的所有条件概率满足

则称该序列具有马尔可夫性。

即 $X_{j}$ 的概率仅依赖于其前一个随机变量 $X_{j - 1}$，而不依赖于其他变量。随机变量的顺序非常重要。

为简单起见，我们仅考虑取值有限且离散的随机变量序列。定义 $P_{i j} = P ( X_{i} | X_{j} )$ 和 $p_{i} = P ( X_{i} )$。马尔可夫链由一个初始概率分布 $p_{i}^{0}$ 和一个转移概率矩阵 $P_{i j}$ 定义，该矩阵满足对于每个 $j$ 有归一化条件 $\textstyle \sum_{i} P_{i j} = 1$。通常链通过递归方式生成：

其中 $p_{i}^{( 0 )} = p_{i}^{0}$。如果对等式两边关于 $i$ 求和，并利用归一化条件 $\textstyle \sum_{i} P_{i j} = 1$ 和 $\textstyle \sum_{i} p_{i}^{0} = 1$，我们会发现 $\begin{array} {r} {\sum_{i} p_{i}^{( k )} = 1} \end{array}$，因此迭代过程保持概率守恒。

方程(2.17)的显著之处在于，在关于 $P_{i j}$ 的非常一般的条件下，且与起始点 $p_{i}^{0}$ 无关，迭代最终会达到一个满足平稳条件的状态 $p_{i}$：

也就是说，它产生了矩阵 $P_{i j}$ 的特征向量 $p_{i}$，其特征值为 1，且向量分量满足 $\textstyle \sum_{i} p_{i} = 1$。

一般而言，$p_{i}$ 是未知的。正如我们将看到的，在基态量子蒙特卡罗模拟中，找到 $p_{i}$ 是目标。在其他应用中，如伊辛模型模拟，则必须生成特定的分布。为此需要选择合适的转移概率。Metropolis 等人 (1953) 提出，如果转移概率满足下式，我们将得到特定的分布 $p_{i}$：

这就是细致平衡条件，它使得处于状态 $j$ 并转移到 $i$ 的概率与处于状态 $i$ 并转移到 $j$ 的概率相等$^{3}$。如果将(2.19)的右边代入(2.18)的右边，则 $\begin{array} {r} {p_{i} = \sum_{j} P_{j i} p_{i} = p_{i}} \end{array}$，这里我们使用了归一化条件 $\textstyle \sum_{j} P_{j i} = 1$。因此，在平稳状态下，$P_{i j}$ 和 $p_{i}$ 是一致的。正如我们即将讨论的，Metropolis 等人还给出了这种 $P_{i j}$ 的一般定义以及从中采样的算法。该算法称为 Metropolis 算法，它具有一个显著特性：被采样的函数不需要归一化。

2.4.2 随机矩阵（Stochastic Matrices）#

现在，我们准备定义蒙特卡罗算法，这些算法使我们能够回答关于相互作用粒子系统的两个基本问题，即：其基态性质是什么？以及其有限温度性质是什么？为了找到基态，我们需要构造一个马尔可夫链，其转移概率能够投影到基态上去。这正是费米所建议的那种蒙特卡罗计算（第1章）。为了计算热力学量，我们需要一个满足玻尔兹曼分布细致平衡条件的转移概率。我们将基态算法的讨论推迟到第9、10和11章。在我们开始讨论用于平衡态热力学的几种细致平衡算法之前，首先证明马尔可夫链的非凡性质，并讨论其有效性的普遍条件。

马尔可夫链总是存在一个满足(2.18)的唯一平稳分布，这一性质源于转移概率 $P_{i j}$ 是一个遍历的随机矩阵（ergodic stochastic matrix）。随机矩阵是指元素非负（非负矩阵）且每列之和等于 1 的矩阵，即 $\textstyle \sum_{i} P_{i j} = 1^{4}$。不可约性意味着，不可能通过一系列行和列的置换将 $P_{i j}$ 变换成以下形式：

其中 $A$ 和 $C$ 是任意阶数的方阵。如果存在这样的置换，则称该矩阵为可约矩阵（reducible）。类似地，非周期性（aperiodicity）指不存在一个置换可将 $P_{i j}$ 转化为如下形式：

一个既不可约又非周期的随机矩阵（stochastic matrix）称为遍历矩阵（ergodic matrix）。⁵

一般而言，随机矩阵不是对称的，即 $P_{i j} \neq P_{j i}$ 。因此，它具有共享相同特征值的不同左特征向量 $( x^{\alpha} )$ 和右特征向量 $( \boldsymbol{y}^{\alpha} )$：

如果特征值互异，$\lambda_{\alpha} \neq \lambda_{\beta}$ ，则右特征向量 $x^{\alpha}$ 和左特征向量 $y^{\beta}$ 线性无关，且满足 $[ x^{\alpha} ]^{T} \cdot y^{\beta} = 0$。

虽然证明超出了本书范围，但可以证明，一个不可约随机矩阵有一个等于1的非退化特征值（nondegenerate eigenvalue），且对应的右特征向量的所有分量为正（Meyer, 2000, 第8章）。这些结果来自Perron-Frobenius定理的应用。⁶ 该定理指出，如果一个不可约矩阵（不一定是随机的）是非负的，那么它有一个等于其谱半径（spectral radius）的非退化特征值。⁷ 因此，该特征值是实数且为正。定理进一步指出，相应的右特征向量的分量也为正，并且该特征向量是唯一的（忽略整体缩放因子）。如果矩阵是随机的，那么我们可以证明特征值的绝对值不能超过1，因此至少有一个占优特征值（dominating eigenvalue）为1，其他占优特征值（如果有的话）是绝对值为1的复数。

为了进一步确定占优特征值，我们需要附加一个条件，即非周期性。我们可以证明，如果矩阵 $P$ 既是不可约的又是非周期的，那么存在一个有限数 $n$ 使得 $P^{n}$ 是正矩阵（positive matrix）；也就是说，对于任意 $i$ 和 $j$，有 $[ P^{n} ]_{i j} > 0$。当一个非负矩阵满足此条件时，称其为原始矩阵（primitive），相应的采样称为遍历（ergodic）。

如果一个随机矩阵是正矩阵，那么从任何给定状态出发，链可以在一步内转移到任何其他状态。注意，一个不可约随机矩阵通常含有零元素（并且可能有多个模为1的特征值）。这意味着在马尔可夫链中，某些状态无法从给定状态一步到达。⁸ 进一步可以证明（Meyer, 2000, 第8章）：对于正随机矩阵，只有一个特征值为1的特征向量。很容易证明这些性质同样适用于原始随机矩阵。也就是说，如果一个随机矩阵是原始的，则 $P$ 的特征值为1的右特征向量唯一且为正。此外，当 $P$ 的选择满足给定目标分布 $p$ 的平稳条件（2.18）时，右占优特征向量等于 $p$。这意味着在马尔可夫链蒙特卡洛模拟中经过足够多次迭代后，概率分布收敛到目标分布。

现在，取某个（未归一化的）向量 $\psi$，将其表示为右特征向量 $y^{\alpha}$ 的线性组合，即 $\begin{array} {r} {\psi = \sum_{\alpha} a_{\alpha} y^{\alpha}} \end{array}$ ，其中 $a_{1} = 1$ ，并按 $1 > | \lambda_{2} | \geq | \lambda_{3} | \geq ...$ 对特征值排序。将该向量与矩阵 $P$ 反复相乘得到

当$k$变大时，具有次优特征值的特征向量被投影掉，只剩下非简并单位特征值的右特征向量。令$y_{i} = p_{i} > 0$便得到(2.18)。我们注意到这个投影状态不依赖于起始状态$p_{i}^{0}$。这个练习揭示了马尔可夫链蒙特卡洛（Markov chain Monte Carlo）的另一个重要特性：在开始采样之前，我们必须迭代一段时间以“平衡”系统。这与之前讨论的直接采样方法形成对比。

虽然一个算法在形式上可能是遍历的（ergodic），但在实际模拟中它可能表现得仿佛并非如此。例如，相空间中可能存在几个高概率区域，它们由低概率路径连接。在这些高概率区域之间移动是一种罕见事件，会耗费大量计算时间。结果，模拟可能锁定在其中一个区域，导致虚假收敛，并产生不能代表整个相空间的构型。在这种情况下，正确的采样需要很长的模拟时间，并且数据分析（第3章）需要特别小心以确保测量值之间的统计独立性。

2.5 细致平衡算法#

我们现在讨论两类满足细致平衡（detailed balance）的算法，因此允许从预先指定的概率分布中进行采样。第一类包括Metropolis算法及相关算法；第二类是热浴算法（heat-bath algorithm）及相关算法。对于这两类算法，新状态通常通过相对局部的构型变化来采样。与Metropolis算法相比，热浴算法没有明确的拒绝步骤。在更广泛的文献中，该算法被称为吉布斯采样器（Gibbs sampler），特别是在统计学文献中，存在大量吉布斯采样器的变体。尽管这两类算法采用不同的方法论，但热浴算法实际上是所谓Metropolis-Hastings算法的一个特例。

2.5.1 Metropolis算法#

Metropolis算法（Metropolis等人，1953）非常巧妙，其优雅源自其简单性。它最初被提出时，更多被视为一种“看起来合理”的方法，而非与马尔可夫链数学有明确联系的方法。发明者确实认识到它是一种计算工具，原则上可以解决经典平衡统计力学中的任何问题。我们引用Metropolis等人（1953）引言段落的原文：

本文的目的是描述一种通用方法，适用于快速电子计算机，用于计算任何由相互作用的单个分子组成的物质的性质。假设经典统计成立……

事实上，Metropolis算法已经持续了几十年，其应用远远超出了统计力学领域，扩展到所有使用随机过程的分析领域。我们现在介绍该算法，并证明它定义了一个满足细致平衡条件的随机转移矩阵。然后我们解释Metropolis等人从该矩阵采样的过程，并讨论其在统计力学中的应用。

为了渐近地采样概率分布$p_{i}$，Metropolis等人提出了以下形式的转移概率矩阵

其中$T_{i j}$是试转移（提议）概率的对称矩阵的元素，满足

而$A_{i j}$是接受矩阵的元素，

通常写作

严格来说，定义(2.22)和(2.23)仅对$i \neq j$成立。这一条件通常是算法实现中根深蒂固的部分：通过使用$T_{ij}$提出的内容本质上是一种变化，因此$i \neq j$。然而，为了证明该算法定义的转移概率矩阵是随机矩阵，让我们将转移概率矩阵表述为

$P_{jj}$的表达式定义了Metropolis算法的拒绝分量。如果提出的转移被拒绝，则当前状态会被添加到马尔可夫链中。

现在转移概率矩阵已完全定义，其随机性很容易证明。⁹ 我们只需证明$\textstyle \sum_{i} P_{i j} = 1$，因为从其定义中明显可知$P_{i j} \ge 0$。证明步骤如下：

为了证明详细平衡条件(2.19)，我们首先注意到它对$i = j$的情况是显然成立的。然后，对于$i \neq j$，在不失一般性的前提下假设$p_i < p_j$（由于标签的任意性），

由此可得 $P_{i j} p_{j} = P_{j i} p_{i}$ 。在得出这一结论时，我们利用了T的对称性、(2.20)和(2.23)。至此，详细平衡（detailed balance）对所有i和j成立这一证明完成。

算法3：Metropolis算法 输入：提议概率T，极限分布 $p$，配置j。 从 $T_{i j}$ 中采样一个i 生成一个均匀随机数 $\zeta \in [0,1]$ 如果 $\zeta \leq p_i / p_j$，则 $j \leftarrow i$； 结束条件 返回 $j$

给定j时采样i的过程异常简单（见算法3），但事后看来可能并不显而易见。在此过程中，Metropolis等人使用 $T_{i j}$ 为马尔可夫链中的下一个配置i生成提议。这个新配置原则上可以是任何允许的配置。一旦提出，他们便利用接受矩阵的 $A_{i j}$ 元素来决定接受i还是j。

这种采样过程在实际中意味着什么？由(2.23)可知，转移概率仅取决于提议配置与当前配置权重之比，因此与其归一化常数无关。在经典统计力学中，此比值为玻尔兹曼因子之比，等于 $\exp(-\Delta E/kT)$ ，其中 $\Delta E$ 是两个配置之间的能量差。马尔可夫链的构建方式为：确定性地或随机性地访问每个格点或粒子，并通过改变该格点上随机变量的状态或移动粒子来提出新配置。这些变化通常是局部的。Metropolis蒙特卡洛模拟很少一次性改变一个配置中的所有变量，甚至不止几个变量。该算法规定：如果新配置具有较低的能量，则“接受提议的局部变化”。此外，该算法还规定“有时接受增加能量的变化，但仅以 $\exp(-\Delta E/kT)$ 的概率接受”。

对于Ising模型，典型的 $T_{i j}$ 会选择特定格点并提议改变自旋状态。该格位上的随机变量，即Ising自旋，有两个可能值：当前值和翻转值。Metropolis转移概率选择这两个配置中的哪一个被添加到马尔可夫链中。值得注意的是，Metropolis算法明确且必然地允许在马尔可夫链中重复同一配置一次或多次。

Metropolis算法是遍历的（ergodic）吗？Metropolis算法将遍历性的负担转嫁给了 $T_{i j}$ 。如果在有限步数内 $T_{i j}$ 能够到达整个相空间，那么Metropolis算法就是遍历的。对于Ising例子，$T_{i j}$ 允许访问每个格点，并且在每次访问时可选择随机变量的任何值。因此，所有自旋配置都可以实现。

2.5.2 广义Metropolis算法#

存在Metropolis算法的推广形式。让我们从一般规则的角度来看待它们：

这里我们不再要求 $T_{i j} = T_{j i}$ 。首先注意到存在一类与接受率相关的通用算法：

其中 $S_{i j}$ 是任意非负对称矩阵¹⁰，它确保对所有i和j都有 $0 \leq A_{i j} \leq 1$ 。$S_{i j}$ 有两种常见选择。其一是

这导致

这被称为Metropolis-Barker算法（Barker, 1965）。第二种选择由Hastings (1970) 提出：

这导致基于以下规则的Metropolis-Hastings算法：

如果 (T) 是对称的，Metropolis-Hastings算法就简化为原始的Metropolis算法。我们注意到，Barker算法和Hastings算法与原始Metropolis算法一样，都不需要平稳分布的归一化。我们将证明这三种算法都满足细致平衡条件的工作留给了练习。

并非所有满足细致平衡的算法都属于 (2.24) 的Metropolis形式。另一类通用算法的接受矩阵形式为 $A_{i j} = S_{i j} / p_{j} T_{i j}$ 。其中 $S_{i j}$ 是任意非负对称矩阵，确保对所有 $i$ 和 $j$ 有 $0 \leq A_{i j} \leq 1$。

在这种情况下，马尔可夫链的转移概率变得与提议概率无关，但依赖于平稳分布的归一化。在原始 Metropolis 等人的见解之上，还存在许多更一般的推广。挑战不在于创造这些推广，而在于决定哪一种最高效。幸运的是，有几个值得注意的指导原则。

首先，我们评论一下计算效率（computational efficiency）与统计效率（statistical efficiency）之间的区别。计算效率通过执行算法每步所需的平均计算机时间来衡量。这里，我们讨论统计效率。从第2.4.2节的讨论中，我们看到，在马尔可夫链收敛的意义上，统计效率由转移概率矩阵的第二大特征值（the second-largest eigenvalue）的大小控制：特征值越小越好。由于矩阵的规模，¹¹ 通常无法定量得知这个特征值。最近，人们开发了蒙特卡洛方法来确定这个特征值 (Rubenstein et al., 2010)。已有的结果支持了 Peskun (1973) 的分析工作。

Peskun 以及后来的其他人 (Liu, 2001) 的数学分析建立了一个定理：如果一个转移矩阵的非对角元素（off-diagonal elements）大于或等于另一个矩阵的对应元素，那么由非对角占优的转移矩阵生成的马尔可夫链所产生的任何期望值的方差，都小于由另一个链产生的方差。Peskun 提出的通用法则是，统计效率随着 $P_{i j}$ 的非对角矩阵元素主导对角元素而增加。这是合理的。非对角占优降低了状态在链中被重复的概率，从而减少了可观测量（observable）的连续估计之间的统计相关性。对于离散空间，Peskun 还证明了 Metropolis-Hastings 算法在统计上比 Metropolis-Barker 算法更高效。

总的来说，细致平衡条件和 Metropolis-Hastings 算法为设计从指定概率密度（probability density）采样的算法提供了灵活的工具。该算法有多种可能的选择，并且可能更优化地控制概率流。通过设计非 Metropolis（或 heat-bath）类型的算法，计算效率得到了显著的提高。该算法的拒绝特性（即重复当前配置一次或多次）是统计效率低下的主要原因，因为它阻碍了生成统计独立的测量值。然而，摆脱 Metropolis 算法并不容易。成功的例子是团簇算法（cluster algorithm）和环回算法（loop algorithm），这些将在第4章、第5章和第6章中讨论。

这些算法的重要创新在于从局部配置变化到全局配置变化的一步。

2.5.3 热浴算法（Heat-bath algorithm）#

热浴算法 (Creutz, 1980) 涉及局部移动，如同 Metropolis 算法一样，但与 Metropolis 算法不同的是，对局部变量 $x_{i}$ 值的提议改变独立于其当前值，并且总是被接受。它可能会在马尔可夫序列中被重复多次，但这不是拒绝的结果——而是选择的结果。

策略是从以少量变量为条件的单变量密度中采样 $x_{i}$。换句话说，我们从某个 $P ( x_{i} | x_{i_{1}} , \ldots , x_{i - 1} , x_{i + 1} , \ldots , x_{i_{m}} )$ 中采样 $x_{i}$，其中 $m \ll N$。如果 $x_{i}$ 的局部环境被替换为一个由条件变量定义的热浴（heat bath），那么我们将从这个环境预期的平衡分布中采样。

我们可以通过使用标准提议概率并反复更新同一位置的构型，在Metropolis算法中实现相同的效果。最终，我们将开始从局部平衡概率中采样。因此，热浴算法（heat-bath algorithm）能在一步之内完成Metropolis算法可能需要很多步才能完成的事情。正因如此，热浴算法通常被认为比Metropolis算法更高效。然而，由于选择导致的链中状态重复，可能会产生与Metropolis情况下拒绝相同的效果。关于哪种算法更高效的情况实际上更为复杂（Rubenstein等，2010）。

为了说明该算法，我们再次以Ising模型为例。我们选择某个格点 $i$，然后固定其他格点上自旋的值。所选格点自旋的局部能量为 $\begin{array} {r} {E ( s_{i} ) = - J s_{i} \sum_{j \neq i} s_{j}} \end{array}$，从而得到局部玻尔兹曼因子（Boltzmann factor）$\exp [ - E ( s_{i} ) / k T ]$。考虑到两种自旋取向，我们可以将自旋 $s_{i} = \pm 1$ 的局部概率写为

从这个概率中采样，会在位置 $i$ 处放置一个值为 $+ 1$ 或 $- 1$ 的 $s_{i}$，从而从旧构型生成一个新构型，而不考虑所选格点上自旋的先前值。

为了提供更多细节和理由，我们按如下步骤进行。如果一组随机变量的概率函数为 $P ( s_{1} , s_{2} , \ldots , s_{N} )$，热浴算法通过固定其他变量，采样每个变量的条件概率，从旧构型生成新构型（见算法4）。

算法4 热浴算法（Heat-bath algorithm）。 输入：一个构型 $( s_{1} , s_{2} , \ldots , s_{N} )$ 对于 $i = 1$ 到 $N$ 执行： 从 $P ( s_{i}^{\prime} | s_{1} , \ldots , s_{i - 1} , s_{i + 1} , \ldots , s_{N} )$ 中采样 $s_{i}^{\prime}$ $s_{i} \gets s_{i}^{\prime}$； 循环结束 返回更新后的 $( s_{1} , s_{2} , \ldots , s_{N} )$

多元条件概率的定义

加上边缘化（marginalization）的定义

给出了我们所需的概率。对于具有最近邻相互作用的零场一维Ising模型，我们可以用一个简单公式来表达这些概率。条件概率变为

这是式（2.30）的一个更详细表述。注意，我们从所有的随机变量出发，但最终得到的表达式只依赖于少数几个变量。这种简化是伊辛模型短程交换相互作用的结果。

热浴算法是 Metropolis-Hastings 算法（2.28）的一个特例，其接受概率等于 1。为了看清这一点，让我们在热浴算法的第 $i$ 步构造 Metropolis-Hastings 接受比

这里我们两次使用了（2.31）式，一次用于 $P ( s_{1} , \ldots , s_{i}^{\prime} , \ldots s_{N} )$，另一次用于 $P ( s_{1} , \ldots , s_{i} , \ldots s_{N} )$。因此，热浴算法的每一步都对应于一个满足 $A_{i j} = 1$ 的 Metropolis-Hastings 抽样。

2.6 Rosenbluth 定理（Rosenbluth's theorem）#

如前所述，细致平衡和 Metropolis 算法最初更多是被当作看似合理的东西提出，而非作为数学推导的结果。在 Metropolis 等人发表论文大约四年后，Wood 和 Parker（1957）将该算法与马尔可夫过程联系起来。这一虽然经常被引用但很少被分析的论述，为原始提案提供了标准的合理性论证。在 Metropolis 等人的论文之后不久，Rosenbluth（1953）撰写了一份未发表、且直到最近才被注意到的报告（Gubernatis, 2003），证明了该算法的有效性。他的证明聚焦于概率的守恒。事实上，他证明了在遍历性和细致平衡的条件下，概率流在相空间中的流动方式使得概率密度与正则系综之间的平均平方偏差不仅变为零，而且是单调地变为零。他的深刻证明指出了收敛过程的一个特殊性质，这是标准马尔可夫链证明中所没有的。

这里我们提供一个相关但更简单的证明，用于证明马尔可夫链的收敛性和 H 定理（H-theorem）。虽然在一般背景下证明它们需要 Perron-Frobenius 定理（Perron-Frobenius theorem），而该定理本身就需要很多行来证明，但假设平稳条件（2.18）被目标分布满足，则简化了任务。大多数蒙特卡洛算法满足一个更强的条件，即细致平衡条件，这一事实证实了该假设的合理性。如我们将看到的，在证明马尔可夫过程的收敛性时，我们只需要平稳条件和遍历性（第 2.4.2 节）。我们将在建立 H 定理或广义 Rosenbluth 定理（generalized Rosenbluth's theorem）时需要细致平衡条件（2.19）。

我们考虑两个系综 $\mathcal{E}_{p}$ 和 $\mathcal{E}_{q}$，它们的状态 $i$ 的概率密度分别为 $p_{i}$ 和 $q_{i}$。我们定义这些系综之间的距离为

现在考虑 $p_{i}$ 和 $p_{i}^{\prime}$ 通过马尔可夫链的一步相关联的情况：$p_{i}^{\prime} = \sum_{j} P_{i j} p_{j}$。我们已经注意到，如果 $P$ 是一个随机矩阵，那么这个操作能保持概率守恒。系综 $\mathcal{E}^{\prime}$ 与满足平稳条件的目标系综 $\mathcal{E}^{\mathrm{eq}}$ 之间的距离为

注意我们在第一步中使用了平稳条件 $\begin{array} {r} {\sum_{j} P_{i j} p_{j}^{\mathrm{eq}} = p_{i}^{\mathrm{eq}}} \end{array}$，并在最后一步中使用了概率守恒条件 $\textstyle \sum_{i} P_{i j} = 1$。¹² 在右侧，我们识别出 $\| \mathcal{E} - \mathcal{E}^{\mathrm{eq}} \|$，因此

因此，算法的每一步都减小了当前系综与平衡系综之间的距离。由于该距离显然以零为下界，它必定收敛到某个值。剩下的问题是这个值是否为零。我们可以通过更仔细地检查 (2.32) 来回答这个问题。

首先，我们注意到，即使将 (2.32) 中的 $P_{i j}$ 替换为 $[ P^{n} ]_{i j}$，并让 $\mathcal{E}^{\prime}$ 为 n 步之后的系综，(2.32) 也必然成立。令 n 足够大，使得 $[ P^{n} ]_{i j} > 0$ 对所有 i 和 j 的组合成立。由于遍历性（ergodicity），我们可以选择这样的 n。现在，当收敛达到后，(2.32) 中的等式必须成立。为了使等式成立，对于所有从同一 i 出发经 n 步可到达的 j，$p_{j} - p_{j}^{\mathrm{eq}}$ 必须具有相同的符号（或为零）。由于对所有 j 都有 $[ P^{n} ]_{i j} > 0$，这意味着要么对所有 j 有 $p_{j} \leq p_{j}^{\mathrm{eq}}$，要么对所有 j 有 $p_{j} \geq p_{j}^{\mathrm{eq}}$。但由于归一化条件 $\begin{array} {r} {\sum_{j} p_{j} = \sum_{j} p_{j}^{\mathrm{eq}^{\cdot}} = 1} \end{array}$，每个不等式都必须化简为 $p_{j} = p_{j}^{\mathrm{eq}}$。因此，$p_{j}$ 必须收敛到 $p_{j}^{\mathrm{{e q}}}$。这一观察完成了定理的证明，即具有平稳条件（或细致平衡条件）的遍历马尔可夫过程必然收敛到目标分布。由此也可推出，模最大的特征值为 1，且对应的特征向量唯一且等于 $p_{i}^{\mathrm{eq}}$，否则该过程将不会总是从任意初始分布收敛到 $p_{i}^{\mathrm{eq}}$。

我们还可以证明，如果马尔可夫过程满足细致平衡，则自由能单调递减并最终收敛到其平衡值。我们称这种蒙特卡罗模拟中的 H 定理为 Rosenbluth 定理（Rosenbluth's theorem）。这种单调性是比平衡密度的平均平方偏差更一般的凸函数的特征性质（Renyi, 1960; Kawashima, 2007）。断言如下：

结合条件 $P_{i j} p_{j}^{\mathrm{eq}} = P_{j i} p_{i}^{\mathrm{eq}}$，导致

单调收敛，如果 $f ( x )$ 是一个凸函数。这一类可观测量最重要的例子可以通过令 $f ( x ) = x \ln x$ 得到。结果是 Kullback-Leibler 信息（Kullback-Leibler information）：

对于物理学家来说，这个量更熟悉的解释是过剩自由能（excess free energy），它度量当前自由能比其平衡值大多少。即，当 $p_{i}^{\mathrm{eq}} \propto \stackrel{\mathrm{-}} {e}^{- E_{i} / k T}$ 时，$I_{\mathrm{KL}} = ( F - F^{\mathrm{eq}} ) / ( k_{B} T )$，其中 $F = \langle E \rangle - T S$ 和 $F^{\mathrm{eq}} = \langle E \rangle^{\mathrm{eq}} - T S^{\mathrm{eq}}$，S 是冯·诺依曼熵（von Neuman entropy），$S \equiv - \sum p \log p$。如果我们令 $f ( x ) = x^{n}$，我们得到 n 阶 Renyi 信息（Renyi information）：

证明如下。在迭代的某一步，我们有

利用 (2.33)，我们将该方程重写为

利用细致平衡条件（detailed balance condition），该式变为

现在，我们对凸函数 $f$ 运用 Jensen 不等式，即对于任意概率分布 $a_{i}$ 有 $f ( \sum_{i} a_{i} x_{i} ) \leq \textstyle \sum_{i} a_{i} f ( x_{i} )$，从而得到

利用细致平衡条件和 $\textstyle \sum_{i} P_{i j} = 1$，我们发现

这个简洁推导给出了一个与 Rosenbluth 提出的证明截然不同的方法。他的证明是在连续空间中进行的，其论证基础是粒子在相空间中运动时，必须保持微分相空间体积元 dr 中的粒子数 $\rho ( r ) d r$ 守恒。

2.7 熵含量（Entropy content）#

分布的熵，

提供了关于其全局特征的信息。例如，如果所有 $p_{i}$ 都相等，那么对于 n 状态分布，熵取最大值 ln n，我们说该分布包含的信息量最小。换句话说，一个事件与其他任何事件的可能性相同，对于链接下来会产生哪个事件，我们拥有最大的不确定性。如果某个特定事件的概率为 1，这意味着所有其他事件的概率必须为零，那么熵取最小值零。该事件是确定无疑的，没有任何不确定性。一旦马尔可夫链到达平稳（stationarity），即采样一个固定分布，其由所生成分布定义的信息熵含量（information entropy content）就变为常数。我们称这个常数为马尔可夫过程（Markov process）的熵含量。（事实上，在平稳点实际被最小化的信息是 Kullback-Leibler 信息，它正比于过剩自由能（excess free energy）；参见第 2.6 节。）

如果链在到达某平稳点之后又走了 r 步，它就会以概率 $P ( K ) = P_{k_{r} , k_{r - 1}} \cdot \cdot \cdot P_{k_{3} , k_{2}} P_{k_{2} , k_{1}} p_{k_{1}}$ 生成一个事件序列 $K = \{k_{1} , k_{2} , \ldots , k_{r} \}$。让我们考虑所有 $n^{r}$ 个序列的集合。那么，可以证明（Khinchin, 1957）：给定任意小的 $\varepsilon > 0$ 和 $\eta > 0$，对于足够大的 r，该集合被分成两个子集，其中一子集的性质是其内任何序列的概率 P(K) 满足

而另一子集的性质是其内所有序列的概率之和小于 ε。表达式

衡量的是马尔可夫链每前进一步时所获得的信息量。

我们注意到，总共有 $n^{r}$ 种可能的 r 项序列。如果我们将这些序列按概率 $P ( K )$ 递减的顺序排列，那么另一个定理（Khinchin, 1957）指出：如果我们将这些概率相加，直到总和刚好超过满足 $0 < \lambda < 1$ 的正数 λ，那么所用到的序列数目 $N_{r} ( \lambda )$ 满足

与 $\lambda$ 的取值无关

虽然这两个定理的证明超出了本书的范围，但它们指出了重要的一点：当可能的事件序列数量等于 $n^{r} = \exp ( r \ln n )$ 时，马尔可夫链选择的数量约为 $\exp ( r R )$。换句话说，$e^{R}$ 是我们在马尔可夫链中为下一步选择的有效候选状态数。由于熵的最大值是 ln n，我们几乎可以肯定 $R < \ln n$。例如，在伊辛模型（Ising model）的单自旋更新中，我们随机选择 N 个自旋中的一个，并在两个可能的选择（自旋向上或自旋向下）中为下一个蒙特卡罗步骤（Monte Carlo step）的自旋值做出选择。由于下一个自旋构型有 $N \times 2$ 种选择，R 的上限为 ln(2N)，并且很自然地假设 R 的量级为 ln N。在这个例子中，考虑到 $n = 2^{N}$，R 不仅小于 ln $n$，对于大系统，它远小于 ln n。从这个例子中可以清楚地看到，在典型情况下，辛钦的后一个定理表明，总序列中微不足道的一小部分占据了几乎所有的概率；也就是说，对于这一小部分 $\textstyle \sum_{K} P ( K ) \approx 1$。这就是马尔可夫链采样的力量所在。

我们还可以证明，随着链步离开初始状态，熵会减小。熵的行为虽然相似，但确实与链向平稳状态的实际松弛行为不同，后者取决于转移概率矩阵的第二大特征值的大小。链的熵含量确实为马尔可夫链的行为及其为何能成为大量随机变量分布的有效采样器提供了一个互补的视角。

推荐阅读#

M. H. Kalos and P. A. Whitlock, Monte Carlo Methods, vol. 1: The Basics (New York: Wiley-Interscience, 1986), chapters 1–3 and appendix.

W. H. Press, S. A. Teukolsky, W. T. Vettering, and B. P. Flannery, Numerical Recipes (Cambridge University, 1992), chapter 7.

J. S. Liu, Monte Carlo Strategies in Scientific Computing (New York: Springer, 2001), chapters 5, 7, and 13.

练习#

2.1 使用数学归纳法证明 (2.13) 的有效性；即先证明结果对 $n = 1$ 成立，假设它对 n 成立，然后证明它对 $n + 1$ 也成立。

2.2 如果 $0 \le p \le 0$ 且 $q = 1 - p$，二项分布（binomial distribution）

描述了 n 次试验中 k 次成功的概率。证明采样 k 的一种方法是

2.3 如果 $0 \le p \le 0$ 且 $q = 1 - p$，证明

其中

现在令 $p = \lambda / n$，其中 $n > \lambda$，并证明当 $n \to \infty$ 时，$Q_{n} \to 1$ 且

这个结果表明，对于大 n 和小 $p$，只要 $n p$ 为常数，二项分布近似等于泊松分布（Poisson distribution）。

2.4 提出一种从密度

中采样 x 的方法，分别考虑 $f_{i}$ 易于归一化和不易归一化两种情况。

2.5 考虑一个由随机矩阵

定义的两状态马尔可夫链。

1. 求其特征值，并验证主特征值 $\lambda_{1} = 1$。

2. 求每个 $\lambda_{i}$ 的左、右特征向量 $x_{i}$ 和 $y_{i}$，并验证主特征值的左特征向量分量全部相等，且当 $i \neq j$ 时 $x_{i}^{T} y_{j} = 0$。

3. 利用这些特征向量解析地求出 $P^{n}$。

4. 证明：当 $n \to \infty$ 时，若向量 $p_{0}$ 的分量之和为 1，则 $P^{n} p^{0}$ 与 $p_{0}$ 无关。

5. 证明 $P^{n} p^{0}$ 的极限向量是与 $\lambda_{1}$ 关联的右特征向量。

2.6 证明 Metropolis-Barker 算法（2.26）和 Metropolis-Hastings 算法（2.28）满足细致平衡条件。对热浴算法（heat-bath algorithm）做同样证明。

2.7 典型模拟只使用一种采样方法，例如 Metropolis 方法或热浴方法。对于热浴方法，条件概率并不总是像伊辛模型（Ising model）那样容易采样。

1. 讨论在热浴算法内部使用 Metropolis 算法进行采样的算法问题和机会。

2. 提出一个在 Metropolis 算法内部使用热浴算法可能更有利的场景。

2.8 若提议概率 $T_{i j}$ 与当前状态无关且等于 $\pi ( i )$，则 Metropolis-Hastings 算法的转移概率为

其中 $p_{i}$ 是目标分布，$w_{i} = p_{i} / \pi_{i}$。

1. 证明该转移概率满足细致平衡条件。

2. 若状态按 $w_{1} \geq w_{2} \geq \cdot \cdot \cdot \geq w_{n}$ 标记，证明

其中 $\begin{array} {r} {\lambda_{k} = \sum_{k} ( \pi_{i} - p_{i} / w_{k} )} \end{array}$。利用你对转移矩阵第一特征值大小及特征向量分量的知识，证明 $\lambda_{2} = 1 - 1 / w_{1}$。对应的特征向量是什么？

3. 论证 $\lambda_{k}$ 是当下一步停留在当前状态 $k$ 时的拒绝概率。Liu（2001）已证明，对于该转移概率，$\lambda_{k}$ 是 $P$ 的特征值。

2.9 从熵的定义

和联合熵

出发，证明 $S ( A , B ) = S ( A ) + S ( A | B )$。