实际工作中经常遇到分布复杂的高维随机向量抽样问题。 12的重要抽样法可以应付维数不太高的情况， 但是对于维数很高而且分布很复杂（比如，分布密度多峰而且位置不易确定的情况）则难以处理。

MCMC(马氏链蒙特卡洛)是一种对高维随机向量抽样的方法， 此方法模拟一个马氏链， 使马氏链的平稳分布为目标分布， 由此产生大量的近似服从目标分布的样本， 但样本不是相互独立的。 MCMC的目标分布密度函数或概率函数可以只计算到差一个常数倍的值。 MCMC方法适用范围广， 近年来获得了广泛的应用。

先介绍马氏链的概念（参考 钱敏平 et al. (2011) ）。

设\(\{X_t, t=0, 1, \dots \}\)为随机变量序列， 称为一个随机过程。 称\(X_t\)为“系统在时刻\(t\)的状态”。 为讨论简单起见， 设所有\(X_t\)均取值于有限集合\(S = \{ 1,2,\dots, m\}\), 称\(S\)为状态空间。 如果\(\{ X_t \}\)满足 \[\begin{align} & P(X_{t+1} = j | X_0=k_0, \dots X_{t-1}=k_{t-1}, X_t=i) \nonumber\\ =& P(X_{t+1} = j | X_t=i) = p_{ij}, \ t=0,1,\dots, \ k_0, \dots, k_{t-1}, i, j \in S, \tag{19.1} \end{align}\] 则称\(\{ X_t \}\)为马氏链, \(p_{ij}\)为转移概率， 矩阵\(P = (p_{ij})_{m\times m}\)为转移概率矩阵。 显然\(\sum_{j=1}^m p_{ij}=1, \ i=1,2,\dots,m\)。 对马氏链， \(P(X_{t+k}=j | X_t = i) = p_{ij}^{(k)}\)也不依赖于\(t\)， 称为\(k\)步转移概率。 如果对任意\(i, j \in S\), \(i\neq j\)都存在\(k \geq 1\)使得 \(p_{ij}^{(k)} > 0\)则称\(\{X_t \}\)为不可约马氏链。 不可约马氏链的所有状态是互相连通的， 即总能经过若干步后互相转移。 对马氏链\(\{ X_t \}\)的某个状态\(i\), 如果存在\(k \geq 0\)使得 \(p_{ii}^{(k)} > 0\)并且\(p_{ii}^{(k+1)}>0\)， 则称\(i\)是非周期的。 如果一个马氏链所有状态都是非周期的，则该马氏链称为非周期的。 不可约马氏链只要有一个状态是非周期的则所有状态是非周期的。 对只有有限个状态的非周期不可约马氏链有 \[ \begin{aligned} \lim_{n\to\infty} P(X_n = j | X_0=i) = \pi_j, \ i, j=1,2,\dots, m, \end{aligned} \] 其中\(\{\pi_j, j=1,2,\dots,m \}\)为正常数, 满足\(\sum_{j=1}^m \pi_j=1\)， 称为\(\{ X_t \}\)的极限分布。 \(\{ \pi_j \}\)满足方程组 \[\begin{align} \begin{cases} \sum_{i=1}^m \pi_i p_{ij} = \pi_j, \ j=1,2,\dots,m \\ \sum_{j=1}^m \pi_j = 1, \end{cases} \tag{19.2} \end{align}\] 称满足(19.2)的分布\(\{ \pi_j \}\)为平稳分布或不变分布。 对只有有限个状态的非周期不可约马氏链， 极限分布和平稳分布存在且为同一分布。

如果允许状态空间\(S\)为可列个元素， 比如\(S = \{ 0, 1, 2, \ldots \}\), 极限分布的条件需要更多的讨论。 对状态\(i\)， 如从状态\(i\)出发总能再返回状态\(i\)， 则称状态\(i\)是常返的(recurrent)。 状态\(i\)常返可等价地定义为\(P(X_n=i, \text{i.o.} | X_0=i)=1\)， 其中i.o.表示事件发生无穷多次，\(A_n \text{ i.o.}\)即 \(\cap_{k=1}^\infty \cup_{n=k}^\infty A_n\)。

对常返状态\(i\)， 如果从\(i\)出发首次返回\(i\)的时间的期望有限， 称\(i\)是正常返的。 对于不可约马氏链，只要存在正常返状态则所有状态都是正常返， 这时存在唯一的平稳分布\(\pi\)。 非周期、正常返的不可约马氏链存在极限分布， 极限分布就是平稳分布: \[ \lim_{n\to\infty} p_{ij}^{(n)} = \pi_j, \ \forall i, j \in S \]

设正常返的不可约马氏链的平稳分布为\(\pi\), 设\(h(\cdot)\)是状态空间\(S\)上的有界函数， 则 \[\begin{align} P\left( \lim_{n\to\infty} \frac{1}{n} \sum_{k=0}^n h(X_k) = \sum_{x \in S} \pi_x h(x) \right) = 1 \tag{19.3} \end{align}\] 这类似于独立同分布随机变量平均值的强大数律。 当\(Y\)服从平稳分布\(\pi\)时，上式中的极限就等于\(E h(Y)\)。 如果能设计转移概率矩阵满足正常返、不可约性质的马氏链且其平稳分布为\(\pi(\cdot)\), 则从某个初始分布出发按照转移概率模拟一列马氏链， 由(19.3)可以估计关于\(Y \sim \pi(\cdot)\)的随机变量的函数的期望\(Eh(Y)\)。

如何得到有平稳分布的转移概率矩阵？ 一个充分条件是转移概率满足细致平衡条件。 如果存在\(\{ \pi_j, j \in S \}\), \(\pi_j \geq 0\), \(\sum_{j \in S} \pi_j = 1\), 使得 \[\begin{aligned} \pi_i p_{ij} = \pi_j p_{ji}, \ \forall i \neq j, \end{aligned} \] 称这样的马氏链为细致平衡的(detailed balanced)， 这时\(\{ \pi_j \}\)是\(\{ X_t \}\)的平稳分布。 事实上，若\(P(X_t = i) = \pi_i, i \in S\), 则 \[\begin{aligned} P(X_{t+1} = j) = \sum_i \pi_i p_{ij} = \sum_i \pi_j p_{ji} = \pi_j \sum_i p_{ji} = \pi_j, \ \forall j \in S. \end{aligned} \] 只要再验证转移概率矩阵\((p_{ij})\)满足不可约性和正常返性就可以利用(19.3)估计服从平稳分布的随机变量\(Y\)的函数的期望\(Eh(Y)\)了。

马氏链的概念可以推广到\(X_t\)的取值集合\(\mathcal{X}\)为\(\mathbb R^d\)的区域的情形。 如果各\(\{ X_t \}\)的有限维分布是连续型的， 则(19.1)可以改用条件密度表示, 这时的\(\{ X_t \}\)按照随机过程论中的习惯应该称作马氏过程, 但这里还是叫做马氏链。

如果非周期、正常返的不可约马氏链\(\{ X_t \}\)的平稳分布为\(\pi(x), x \in \mathcal{X}\), 则从任意初值出发模拟产生序列\(\{ X_t \}\)， 当\(t\)很大时，\(X_t\)的分布就近似服从\(\pi\)， 抛弃开始的一段后的\(X_t\)序列可以作为分布\(\pi\)的相关的样本， 抛弃的一段序列叫做老化期。 如果只是为了估计\(Y \sim \pi(\cdot)\)的函数的期望， 可以仅要求马氏链为正常返不可约马氏链。 设\(Y \sim \pi(\cdot)\), \(h(y), y \in \mathcal{X}\)是有界函数, 为估计\(\theta = E h(Y)\)， 用\(X_t, t=k+1, \dots, n\)作为\(\pi\)的样本， 用估计量\(\hat\theta = \frac{1}{n-k} \sum_{t=k+1}^n h(X_t)\)来估计\(\theta\)， 则\(\hat\theta\)是\(\theta\)的强相合估计。 老化期长度\(k\)可以从\(\hat\theta\)的变化图形经验地选取。

这样的估计量\(\hat\theta\)是相关样本的平均值， 无法用原来独立样本的公式估计\(\text{Var}(\hat\theta)\)从而得到\(\hat\theta\)的标准误差。 为了估计\(\text{Var}(\hat\theta)\), 可以采用如下的分段平均法。 把样本\(X_{k+1}, \dots, X_n\)分为\(s\)段， 每段\(r\)个(设\(n-k = sr\))。 设第\(j\)段的\(r\)个\(h(X_t)\)的平均值为\(Z_j, j=1,2,\dots,s\), 设\(\{ Z_j, j=1,2,\dots,s \}\)的样本方差为 \(\hat\sigma^2=\frac{1}{s-1} \sum_{j=1}^s (Z_j - \bar Z)^2\)， 因为\(\hat\theta\)等于\(\{Z_j, j=1,2,\dots,s \}\)的样本均值\(\bar Z\)， 当\(r\)足够大时， 可以认为各\(Z_j\)相关性已经很弱， 这时\(\hat\theta\)的方差可以用\(\hat\sigma^2 / s\)估计。 \(r\)的大小依赖于不同时刻的\(X_t\)的相关性强弱， 相关性越强，需要的\(r\)越大。

以上的方法就是MCMC方法（马氏链蒙特卡洛）。 一般地， 对高维或取值空间\(\mathcal{X}\)结构复杂的随机向量\(X\), MCMC方法构造取值于\(\mathcal{X}\)的马氏链， 使其平稳分布为\(X\)的目标分布。 模拟此马氏链， 抛弃开始的部分抽样值， 把剩余部分作为\(X\)的非独立抽样。 非独立抽样的估计效率比独立抽样低。

MCMC方法的关键在于如何从第\(t\)时刻转移到第\(t+1\)时刻。 好的转移算法应该使得马氏链比较快地收敛到平稳分布， 并且不会长时间地停留在取值空间\(\mathcal{X}\)的局部区域内 （在目标分布是多峰分布且峰高度差异较大时容易出现这种问题）。

Metropolis-Hasting方法(MH方法)是一个基本的MCMC算法， 此算法在每一步试探地进行转移(如随机游动)， 如果转移后能提高状态\(x_t\)在目标分布\(\pi\)中的密度值则接受转移结果， 否则以一定的概率决定是转移还是停留不动。

Gibbs抽样是另外一种常用的MCMC方法，此方法轮流沿各个坐标轴方向转移， 且转移概率由当前状态下用其它坐标预测转移方向坐标的条件分布给出。 因为利用了目标分布的条件分布， 所以Gibbs抽样方法的效率比MH方法效率更高。

设随机变量\(X\)分布为\(\pi(x), x \in \mathcal{X}\)。 为论述简单起见仍假设\(\mathcal{X}\)是离散集合。 算法需要一个试转移概率函数\(T(y|x), x, y \in \mathcal{X}\)，满足 \(0 \leq T(y|x) \leq 1\), \(\sum_y T(y|x) = 1\)， 并且 \[\begin{aligned} T(y|x)>0 \Leftrightarrow T(x|y)>0 . \end{aligned}\]

算法首先从\(\mathcal{X}\)中任意取初值\(X^{(0)}\)。 设经过\(t\)步后算法的当前状态为\(X^{(t)}\), 则下一步由试转移分布\(T(y|X^{(t)})\)抽取\(Y\)， 并生成\(U \sim\)U(0,1), 然后按如下规则转移: \[\begin{aligned} X^{(t+1)} = \begin{cases} Y & \text{如果 } U \leq r(X^{(t)}, Y) \\ X^{(t)} & \text{其它.} \end{cases} \end{aligned}\] 其中 \[\begin{align} r(x, y) = \min \left\{ 1, \frac{\pi(y)T(x|y)}{\pi(x)T(y|x)} \right\}. \tag{19.4} \end{align}\]

在MH算法中如果取\(T(y|x) = T(x|y)\)， 则\(r(x,y) = \min\left(1, \frac{\pi(y)}{\pi(x) } \right)\), 相应的算法称为Metropolis抽样法。

如果取\(T(y|x) = g(y)\)(不依赖于\(x\))， 则\(r(x,y) = \min\left(1, \frac{\pi(y)/g(y)}{\pi(x)/ g(x)} \right)\), 相应的算法称为Metropolis独立抽样法, 和重要抽样有相似之处, 试抽样分布\(g(\cdot)\)经常取为相对重尾的分布。

在MH算法中， 目标分布\(\pi(x)\)可以用差一个常数倍的\(\tilde \pi(x) = C \pi(x)\)代替， 这样关于目标分布仅知道差一个常数倍的\(\tilde \pi(x)\)的情形， 也可以使用此算法。

下面说明MH抽样方法的合理性。 我们来验证MH抽样的转移概率\(A(x,y) = P(X^{(t+1)} = y | X^{(t)} = x)\)满足细致平衡条件。 易见 \[\begin{aligned} A(x, y) = \begin{cases} T(y | x) r(x,y), & y \neq x, \\ T(x | x) + \sum_{z \neq x} T(z | x) [1 - r(x,z)], & y = x, \end{cases} \end{aligned}\] 于是当\(x \neq y\)时 \[\begin{aligned} \pi(x) A(x,y) =& \pi(x) T(y | x) \min \left\{ 1, \frac{\pi(y)T(x | y)}{\pi(x)T(y | x)} \right\} = \min \left\{ \pi(x)T(y | x), \pi(y)T(x | y) \right\} , \end{aligned}\] 等式右侧关于\(x, y\)是对称的， 所以等式左侧把\(x, y\)交换后仍相等。 所以，MH构造的马氏链以\(\{ \pi(x) \}\)为平稳分布。 多数情况下MH构造的马氏链也以\(\{ \pi(x) \}\)为极限分布。

(19.4)中的\(r(x,y)\)还可以推广为如下的形式 \[\begin{aligned} \tilde r(x, y) = \frac{\alpha(x, y)}{\pi(x) T(y | x)}, \end{aligned}\] 其中\(\alpha(x,y)\)是任意的满足\(\alpha(x,y)=\alpha(y,x)\) 且使得\(\tilde r(x,y) \leq 1\)的函数。 易见这样的\(\tilde r(x,y)\)仍使得生成的马氏链满足细致平衡条件。

例19.1 \(X\)的取值集合\(\mathcal{X}\)可能是很大的，以至于无法穷举， 目标分布\(\pi(x)\)可能是只能确定到差一个常数倍。

例如，设 \[\begin{aligned} \mathcal{X} =& \{\boldsymbol x = (x_1, x_2, \dots, x_n) \,:\, \\ & \; (x_1, x_2, \dots, x_n) \text{是} (1,2,\dots,n) \text{的排列，使得}\sum_{j=1}^n j x_j > a \}, \end{aligned}\] 其中\(a\)是一个给定的常数。 用\(|\mathcal{X}|\)表示\(\mathcal{X}\)的元素个数， 当\(n\)较大时\(\mathcal{X}\)是\((1,2,\dots,n)\)的所有\(n!\)个排列的一个子集， \(|\mathcal{X}|\)很大， 很难穷举\(\mathcal{X}\)的元素, 从而\(|\mathcal{X}|\)未知。

设\(X\)服从\(\mathcal{X}\)上的均匀分布， 即\(\pi(\boldsymbol x) = C, \ x \in \mathcal{X}\), \(C = 1/|\mathcal{X}|\)但\(C\)未知。 要用MH方法产生\(X\)的抽样序列。

试抽样\(T(\boldsymbol y | \boldsymbol x)\)如果允许转移到所有的\(\boldsymbol y\)是很难执行的， 因为\(\boldsymbol y\)的个数太多了。 我们定义一个\(\boldsymbol x\)的近邻的概念， 仅考虑转移到\(\boldsymbol x\)的近邻。 一种定义是， 如果把\(\boldsymbol x\)的\(n\)个元素中的某两个交换位置后可以得到\(\boldsymbol y \in \mathcal{X}\)， 则\(\boldsymbol y\)称为\(\boldsymbol x\)的一个近邻， 记\(N(\boldsymbol x)\)为\(\boldsymbol x\)的所有近邻的集合， 记\(| N(\boldsymbol x) |\)为\(\boldsymbol x\)的近邻的个数。 当\(n\)很大时，求\(N(\boldsymbol x)\)也需要从\(C_n^2=\frac12 n(n-1)\)个可能的元素中用穷举法选择。 取试转移概率函数为 \[\begin{aligned} T(\boldsymbol y | \boldsymbol x) = \frac{1}{|N(\boldsymbol x)|}, \ \boldsymbol x, \boldsymbol y \in \mathcal{X}, \end{aligned}\] 即从\(\boldsymbol x\)出发，等可能地试转移到\(\boldsymbol x\)的任何一个近邻上。

因为目标分布\(\pi(\boldsymbol x)\)是常数， 所以这时 \[\begin{aligned} r(\boldsymbol x, \boldsymbol y) = \min\left(1, \frac{|N(\boldsymbol x)|}{|N(\boldsymbol y)|} \right), \end{aligned}\] 即从\(\boldsymbol x\)试转移到\(\boldsymbol y\)后， 如果\(\boldsymbol y\)的近邻数不超过\(\boldsymbol x\)的近邻数则确定转移到\(\boldsymbol y\)， 否则，仅按概率\(|N(\boldsymbol x)| / |N(\boldsymbol y)|\)转移到\(\boldsymbol y\)。 这就构成了对\(X\)抽样的MH算法。

Julia实现：

随机游动MH算法是一种Metropolis抽样方法。 随机游动的步幅\(\sigma\)是重要参数，步幅过大导致拒绝率大， 步幅过小使得序列的相关性太强，收敛到平衡态速度太慢。 一个建议选法是试验各种选法， 使得试抽样被接受的概率在0.25到0.35之间。 见 Liu (2001) § 5.4 P. 115.

例19.3 考虑如下的简单气体模型: 在平面区域\(G = [0,A]\times[0,B]\)内有\(K\)个直径为 \(d\)的刚性圆盘。 随机向量\(\boldsymbol X = (x_1, y_1, \ldots, x_K, y_k)\)为这些圆盘的位置坐标。 分布\(\pi(\boldsymbol x)\)是\(G\)内所有允许位置的均匀分布。希望对\(\pi\)抽样。

先找一个初始的允许位置\(\boldsymbol x^{(0)}\)。比如，把圆盘整齐地排列在左上角。

设已得到\(\boldsymbol x^{(t)}\)，随机选取一个圆盘\(i\)，把圆盘\(i\)的位置试移动到 \((x_i', y_i') = (x_i + \delta_i, y_i + \epsilon_i)\)， 其中\(\delta_i, \epsilon_i\)独立同N(0, \(\sigma^2\))分布。 如果得到的位置是允许的则接受结果，否则留在原地不动。

下面给出R实现。 运行参数：