吉布斯抽样

吉布斯抽样（Gibbs sampling）是马尔可夫链蒙特卡洛的一种构造：每次固定其他坐标，从一个坐标的全条件分布重新抽样。单个坐标的条件分布可能很简单，组合起来却能产生具有依赖关系的联合样本。“全条件”指给定所有其余坐标，而非只给定任意挑选的一部分。

两个二元变量X、Y的联合概率取为下表。这是用于演示算法的合成分布。虽然X和Y各自取1的概率均为二分之一，它们同时为1的概率却是0.4，不能各自独立抛一枚公平硬币来生成这一联合分布。

	Y=0	Y=1	行和
X=0	0.4	0.1	0.5
X=1	0.1	0.4	0.5
列和	0.5	0.5	1

固定一列或一行后重新分配

给定Y=0时，只比较表中第一列，故 $Pr(X=1\mid Y=0)=0.1/0.5=0.2$；给定Y=1时则为 $0.4/0.5=0.8$。由于本例对称，Y在给定X时也有同样的两个概率。一次完整扫描按以下顺序进行： $$X^{\prime }\sim \mathrm{Bernoulli}(0.2+0.6Y),Y^{\prime }\sim \mathrm{Bernoulli}(0.2+0.6X^{\prime }).$$ 第二步必须使用刚抽出的X′。这两次抽样使用独立的新随机数，但最终X′、Y′并不独立；条件概率把它们关联起来。

先固定Y所在的列，再除以该列总概率0.5；更新Y时改为使用新X所在的行。

从 $(X,Y)=(0,0)$ 出发，X′=0的概率为0.8，随后Y′=0的概率也为0.8，所以转到00的概率为0.64。转到01的概率为 $0.8\times 0.2=0.16$；转到10为 $0.2\times 0.2=0.04$；转到11为 $0.2\times 0.8=0.16$。

一轮更新保持联合分布

设有限集合上的目标概率为 $\pi (x,y)$，并只在正概率边缘上使用条件分布。若更新前的联合分布已是π，更新X后得到 $(X^{\prime },Y)$ 的概率为 $$\sum _x\pi (x,y)\pi (x^{\prime }\mid y)={\pi }_Y(y)\pi (x^{\prime }\mid y)=\pi (x^{\prime },y).$$ 这说明只重抽X就保持π。同样的计算证明重抽Y也保持π；把两步接起来，完整扫描仍保持π。证明依赖准确的全条件分布，不能把条件均值或最可能值当作随机抽样。

把状态依次排为00、01、10、11，本例完整扫描的转移矩阵是 $$P=\left(\begin{array}{cccc}0.64& 0.16& 0.04& 0.16\\ 0.16& 0.04& 0.16& 0.64\\ 0.64& 0.16& 0.04& 0.16\\ 0.16& 0.04& 0.16& 0.64\end{array}\right).$$ 旧X会在第一步被替换，所以旧Y相同的两行完全一致。将 $\pi =(0.4,0.1,0.1,0.4)$ 左乘该矩阵，得到的第一分量为 $0.5\times 0.64+0.5\times 0.16=0.4$，其余三项依次为0.1、0.1、0.4，直接核验了平稳性。全部矩阵元为正，因而有限链不可约且非周期，从任意初始状态趋向π。

单坐标条件更新是可逆的，固定顺序把两个更新接起来却未必可逆。本例 $\pi (00)P(00,01)=0.064$，反向流量 $\pi (01)P(01,00)=0.016$，并不相等。保持π不要求每对状态流量都相等。若每步以各二分之一的概率只选择X或Y之一更新，其转移核是两个可逆核的等权混合，则仍可逆。

收敛速度可以逐轮计算

记完整扫描结束后的 $r_t=Pr(Y_t=1)$。若旧Y=0，下一轮Y′=1的概率为 $0.8\times 0.2+0.2\times 0.8=0.32$；若旧Y=1，该概率为 $0.2\times 0.2+0.8\times 0.8=0.68$。因此 $$r_{t+1}=0.32+0.36r_t,r_t=\frac{1}{2}(1-0.36^t)$$ （从00起步时 $r_0=0$）。第一轮的联合概率为 $(0.64,0.16,0.04,0.16)$；按此矩阵再走一轮为 $(0.4864,0.1216,0.0784,0.3136)$。逐轮概率与目标的四项都可比较，不能仅检查两个边缘均值是否接近二分之一。

第一轮与第二轮是对转移矩阵作精确计算得到的分布，不是一条随机轨迹的频率。

二元联合分布的Python实现

代码返回n个完整扫描结束后的二元状态。burn以完整扫描为单位，起点为00；每轮恰好使用两个独立的均匀伪随机数。

Python 3

import random

def binary_gibbs(n, burn=0, seed=26):
    if type(n) is not int or n < 1:
        raise ValueError("n must be a positive integer")
    if type(burn) is not int or burn < 0:
        raise ValueError("burn must be a nonnegative integer")
    rng = random.Random(seed)
    x, y = 0, 0
    result = []
    for step in range(n + burn):
        x = int(rng.random() < 0.2 + 0.6 * y)
        y = int(rng.random() < 0.2 + 0.6 * x)
        if step >= burn:
            result.append((x, y))
    return result

draws = binary_gibbs(20000, burn=50)
states = ((0, 0), (0, 1), (1, 0), (1, 1))
print([draws.count(s) / len(draws) for s in states])

每轮时间为常数，总时间 $O(n+\mathrm{b}\mathrm{u}\mathrm{r}\mathrm{n})$，返回列表占 $O(n)$ 空间。两个条件式是本例专用的概率，其他联合分布必须重新推导全条件，不能只更换数据标签。

支撑约束与同步更新的反例

若错误地用旧X、旧Y同时生成新X、新Y，从00出发的下一轮概率会变为 $(0.64,0.16,0.16,0.04)$。更直接地，把目标π作为输入分布，这种同步更新输出的11概率为 $$0.4(0.2{)}^2+0.1(0.8)(0.2)+0.1(0.2)(0.8)+0.4(0.8{)}^2=0.304,$$ 不再等于0.4。一般并行更新需要额外的条件独立结构或专门设计；直接把固定顺序扫描改成同时执行，会改掉算法。

变式。 让目标只在00、11各放二分之一概率，01、10概率均为零。此时给定Y，X只能等于Y；给定X，Y也只能等于X。从00开始会永远留在00，虽然每次条件更新都正确，π也确实平稳。失败处是链不能在两个正概率状态间移动；把高度依赖的变量合成一块更新，才可能跨越这条支撑约束。

连续变量也可作同样的条件更新。例如二元正态分布的全条件仍为正态，但相关性很强时，单坐标步只能在狭窄区域内缓慢移动。条件分布容易抽样，不能单独保证整条链快速混合。

参考资料

• David A. Levin、Yuval Peres：Markov Chains and Mixing Times，第2版，§3.3 Glauber dynamics，条件重抽样的构造。本文使用独立设计的四格分布说明固定扫描。