马尔可夫链蒙特卡洛

马尔可夫链蒙特卡洛（Markov chain Monte Carlo，MCMC）通过构造一条以目标分布为平稳分布的马尔可夫链，用沿途访问的状态近似计算概率和期望。它适用于目标分布难以直接抽样、但能够比较不同状态相对概率的问题。链上相邻样本通常相关；计算误差时，样本总数不能直接当作独立样本数。

设三个方案的相对权重为3、2、1，收益依次为0、3、6，希望按它们的概率计算平均收益。把权重相加便能直接抽样，所以这个小问题本来不需要MCMC；它的作用是让转移、平衡和抽样误差都能逐项算清。在更大的状态空间中，同样的构造可以只计算当前位置和一个候选位置的权重，避开对所有状态求和。

从相对权重构造转移

把状态记为0、1、2，目标概率为 $\pi =(1/2,1/3,1/6)$。收益函数取 $h=(0,3,6)$，其精确期望是 $0\times (1/2)+3\times (1/3)+6\times (1/6)=2$，后面可用链上收益的平均核对。当前位于状态i时，先以各二分之一的概率提议另外两个状态中的一个j，再以 $$\alpha (i,j)=min(1,\frac{w_j}{w_i})$$ 接受候选；没有接受就保留i。这里 $w=(3,2,1)$。从0提议1的概率是 $1/2$，接受概率是 $2/3$，所以真正转移到1的概率为 $1/3$。同理，从0到2的概率为 $1/6$，余下二分之一留在0。完整转移矩阵为 $$P=\left(\begin{array}{ccc}1/2& 1/3& 1/6\\ 1/2& 1/4& 1/4\\ 1/2& 1/2& 0\end{array}\right).$$ 矩阵第i行是当前位置为i时下一步的分布，每行之和为一。拒绝候选也完成了一步，记录中应当再次写入原状态；若只记录成功移动，访问比例通常就变了。

线上标的是平衡时的概率流 π(i)P(i,j)，并非转移概率。自循环未画出；完整转移概率见矩阵。

计算 $\pi P$，第一分量为 $(1/2+1/3+1/6)/2=1/2$，第二分量为 $1/6+1/12+1/12=1/3$，第三分量为 $1/12+1/12=1/6$。因此 $\pi P=\pi$：若某一步已按π分布，下一步仍按π分布。这叫平稳性。

Metropolis–Hastings的平衡证明

在有限状态集合上，设目标 $\pi (i)>0$，提议矩阵q每行和为一。对 $i\ne j$ 且 $q(i,j)>0$，接受率取 $$\alpha (i,j)=min(1,\frac{\pi (j)q(j,i)}{\pi (i)q(i,j)}).$$ $$P(i,j)=q(i,j)\alpha (i,j).$$ 若反向提议概率为零，接受率就是零。未提议的转移概率也为零；对角元定义为 $P(i,i)=1-\sum _{j\ne i}P(i,j)$。把接受率代入，可得 $$\begin{array}{cc}& \pi (i)P(i,j)\\ & =min\{\pi (i)q(i,j),\pi (j)q(j,i)\}\\ & =\pi (j)P(j,i).\end{array}$$ 这个等式称为细致平衡。对i求和，并利用P第j行之和为一，便有 $\sum _i\pi (i)P(i,j)=\pi (j)$，证明目标分布平稳。若只知道未归一化权重 $\pi (i)=w_i/Z$，公共的Z在比值中约去；对称提议的q也相消，回到开头的接受率。

平稳性之外还要检查链能走到哪里。在有限状态情形，若链不可约，即任意两个状态能以正概率经有限步互达，则平稳分布唯一；再加上非周期性，各步分布从任意初始分布趋向π。上例任意两点直接互达，且0可以原地停留，满足这些条件。一般连续或无限状态空间还需相应的遍历性条件，不能照搬有限矩阵的论证。细致平衡是保证平稳的充分条件，固定顺序的吉布斯抽样可以平稳而不满足细致平衡。

单步分布与时间平均

另取一条可以写出精确公式的两状态链： $$P_2=\left(\begin{array}{cc}3/4& 1/4\\ 1/2& 1/2\end{array}\right),{\pi }_2=(2/3,1/3).$$ 从0出发，记 $p_t=Pr(X_t=1)$。按前一步的两个可能状态分解， $$p_{t+1}=(1-p_t)/4+p_t/2=1/4+p_t/4.$$ 减去平衡值 $1/3$，每一步的差缩成原来的四分之一，因此 $$p_t=\frac{1}{3}(1-4^{-t}).$$ 前四个值为 $0,1/4,5/16,21/64$。这些数描述重复运行整条链后，第t步落在1的比例；一次运行只有0或1，不能把一条轨迹的起伏误读成这条平滑收敛规律。

点只对应整数步。青点是精确边际概率，金色虚线是平衡概率；连线用于引导阅读。

实际计算用时间平均 ${\hat{\mu }}_n=n^{-1}\sum _{t=1}^{n}h(X_t)$ 估计 $\mu =\sum _i\pi (i)h(i)$。有限不可约链的遍历定理保证该平均几乎必然趋向μ；这条时间平均结论不要求非周期性。例如必定交替的0、1链，从0出发时单步分布不收敛，但1的时间比例仍趋向二分之一。定理的完整证明见所引马尔可夫链资料；上面的递推已独立证明了本例单步分布的收敛。

若两状态链从平稳分布开始，令 $H_t={\mathrm{𝟏}}_{\{X_t=1\}}$，则 $\mathrm{Var}(H_t)=2/9$，由同一递推可得 $\mathrm{Cov}(H_t,H_{t+k})=(2/9)4^{-k}$。展开平均值的方差，距离k的有序交叉项共有 $2(n-k)$ 个，所以 $$\mathrm{Var}({\hat{\mu }}_n)=\frac{2}{9n^2}[n+2\sum _{k=1}^{n-1}(n-k)4^{-k}].$$ 把它与n个独立伯努利样本的方差 $2/(9n)$ 比较，比例随n趋向 $1+2\sum _{k\ge 1}{4}^{-k}=5/3$。本例的大样本有效样本量约为 $3n/5$。这个值针对所选H和本条链；其他函数、链和非平稳起点都需要另算。

纵轴为相关样本的方差除以同数目独立样本的方差，极限为5/3；数值越大，说明相关性增加的误差越明显。

三状态采样的Python实现

下面只实现开头权重3、2、1和对称提议的例子。函数记录每次转移后的状态，包含拒绝造成的重复；舍去前burn步之后保留n步，不修改外部输入。

Python 3

import random

def sample_three_states(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)
    weights = (3, 2, 1)
    state, samples = 0, []
    for step in range(n + burn):
        candidate = (state + 1 + rng.randrange(2)) % 3
        if rng.random() < min(1.0, weights[candidate] / weights[state]):
            state = candidate
        if step >= burn:
            samples.append(state)
    return samples

draws = sample_three_states(20000, burn=50)
print([draws.count(i) / len(draws) for i in range(3)])
payoff = (0, 3, 6)
print("mean payoff:", sum(payoff[s] for s in draws) / len(draws))

时间代价为 $O(n+\mathrm{b}\mathrm{u}\mathrm{r}\mathrm{n})$，存储返回列表占 $O(n)$；若只累计收益均值，可以只保留计数和总和。固定种子使这一实现可复现，频率应围绕 $(1/2,1/3,1/6)$ 起伏，不能要求一次运行恰好等于它。

丢弃一些初始步有助于减轻起点影响，但没有适用于所有链的固定预热步数；增加预热也不能修复不可达的状态。多起点轨迹、自相关和有效样本量是诊断线索，不能单凭“曲线看起来稳定”证明已充分探索目标分布。

变式。 把权重改为6、4、2，接受率完全不变，目标概率也不变。若把提议改成每次只能沿0→1→2→0的单向环前进，却仍使用一般Metropolis–Hastings公式，则每次候选的反向提议概率为零，所有移动被拒绝。此时转移矩阵为 $P=I$，目标分布虽然平稳，链却可约：从任一起点出发都被困在原状态，不能靠增加运行次数完成所需抽样。

参考资料

• David A. Levin、Yuval Peres：Markov Chains and Mixing Times，第2版，第1章的平稳分布与可逆性、第4章的收敛，以及附录C §C.1、定理C.1的遍历定理（印刷页390）。
• 相关方法：重要性抽样用密度比加权；拒绝抽样通过接受规则得到独立样本。它们与本篇的相关链样本有不同的误差分析。