MCMC之前尘后事

MCMC即马尔可夫链蒙特卡罗（Markov Chain Monte Carlo, MCMC），是一类在统计学和机器学习等领域广泛应用的强大算法。

前尘

最早的源于Metropolis等人发表于《Journal of Chemical Physics》的主要针对以下形式的积分：
$$
J=\int F(\theta)exp(-E(\theta)/kT)d\theta / \int exp(-E(\theta)/kT)d\theta
$$
$\theta$表示$R^2$上一组N个粒子（particles），能量E为：
$$
E(\theta)=\frac{1}{2}\sum _ {i=1}^N\sum _ {j=1,j\ne i}^NV(d _ {ij})
$$
其中V是势函数，d是粒子i和j在$\theta$中的欧几里得距离。玻尔兹曼分布$exp(-E(\theta))$。

考虑到问题的维度很大，即使是标准的蒙特卡洛技术也无法正确地近似$J$。

因为对于粒子系统（在 2N 的正方形中均匀分布）的随机构型的大多数实现，exp{−E(θ)/kT }非常小。为了提高蒙特卡洛方法的效率，Metropolis 等人（1953 年）提出了对 N 个粒子的随机游走修改，即
$$
\begin{align}
x_i’=x_i+\sigma\xi _ {1i}
\\
y_i’=y_i+\sigma\xi _ {2i}
\end{align}
$$
$\xi _ {1i}、\xi _ {2i}$都是均匀分布U(-1,1)。

然后计算新配置（configuration）和旧配置之间的能量差∆E，并以概率接受新配置。
$$
min(1,exp(-\Delta E/KT))
$$
否则，先前配置将被复制，即在随机游走 τ 步的 F (θ) 的最终平均中，其计数器增加一 。但Metropolis 等人一次移动一个粒子，而不是一起移动所有粒子，这使得初始算法看起来像是一种原始的吉布斯采样器。

基本采样方法

像变分贝叶斯推断是基于确定性近似的推断⽅法，考虑基于数值采样的近似推断⽅法，也被称为蒙特卡罗方法。

MCMC方法巧妙地将两个核心概念结合起来：蒙特卡罗方法和马尔可夫链。

• 蒙特卡罗方法：这是一种通过大量随机抽样来估计复杂问题数值解的方法。

• 马尔可夫链：一个马尔可夫链是一个随机过程，其中未来的状态仅依赖于当前状态，而与过去的状态无关（即“无记忆性”）。在MCMC中，我们构造一个特殊的马尔可夫链，使其平稳分布恰好是我们想要抽样的目标概率分布 $\pi(x)$。这意味着，当马尔可夫链运行足够长的时间后，它所产生的样本将近似服从目标分布 $\pi(x)$。

在介绍MCMC之前，我们先介绍一些基本概念和基本方法。

我们希望解决的基本的问题涉及到关于⼀个概率分布p(z)寻找某个函数f(z)的期望。即:
$$
E[f]=\int f(z)p(z)dz
$$
采样⽅法背后的⼀般思想是得到从概率分布p(z)中独⽴抽取的⼀组变量$z^{(l)}$。这使得期望可以通过有限和的⽅式计算，即
$$
\hat f=\frac{1}{L}\sum f(z^{(l)})
$$
且估计的精度不依赖于z的维度，且原则上数量较少也能达到较高的精度。

hoeffding不等式：

一系列随机有界独立变量，$z_1…z_n$，$z_i\in[a,b]$，对所有的i，有：
$$
P(\frac{1}{n}\sum(z_i-E(z_i))\ge t)\le exp(-\frac{2nt^2}{(b-a)^2})
$$

基本采样

我们考虑如何从简单的⾮均匀分布中⽣成随机数，假定我们已经有了⼀个均匀分布的随机数的来源。假设z在区间(0*,* 1)上均匀分布，我们使⽤某个函数f(·)对z的值进⾏变换，即y=f(z)。y上的概率分布为:
$$
p(y)=p(z)\left|\frac{dz}{dy}\right|
$$
我们的⽬标是选择⼀个函数f(z)使得产⽣出的y值具有某种所需的具体的分布形式p(y)，对上式进行积分有：
$$
z=h(y)=\int _ {-\infty}^y p(\hat y)d\hat y
$$
它是p(y)的不定积分。因此，$y=h^-1(z)$，因此我们必须使⽤⼀个函数来对这个均匀分布的随机数进⾏变换

假如指数分布（exponential distribution）$p(y)=\lambda exp(-\lambda y)$，其中$0\le y<\infty$，故有
$$
\begin{align}
z&=h(y)=\int_0^y p(y)dy
\\
&=\int_0^y \lambda exp(-\lambda y)dy
\\
&=1-exp(-\lambda y)
\end{align}
$$

故我们均匀生成z，并通过$y=-\lambda^{-1}ln(1-z)$来生成指数分布。

另外统计学中还有一个Box-Muller⽅法⽤于⽣成⾼斯概率分布的样本。

但是我们这种方法只能用于⼀些⾮常有限的简单的概率分布可⾏，因此我们必须寻找⼀些更加⼀般的⽅法。

拒绝采样

为了应⽤拒绝采样⽅法，我们需要⼀些简单的概率分布q(z)，有时被称为提议分布（proposal distribution）。满足对所有的z，都要$kq(z)\ge \tilde p(z)$，函数kq(z)被称为⽐较函数。

⾸先，我们从概率分布q(z)中⽣成⼀个数z0。

接下来，我们在区间[0, kq(z0)]上的均匀分布中⽣成⼀个数u0。这对随机数在函数kq(z)的曲线下⽅是均匀分布。

最后，如果$u_0>\tilde{p}(z_0)$，那么样本被拒绝，否则u0被保留。因此，如果它位于下图的灰⾊阴影部分，它就会被拒绝。这样，剩余的点对在曲线$\tilde{p}(z)$下⽅是均匀分布的，因此对应的z值服从概率分布p(z)。

在某些情况下，找到满足对所有的z，都要$kq(z)\ge \tilde p(z)$的q是困难的，所以又有可调节的拒绝采样。

拒绝采样在⼀维或⼆维空间中是⼀个有⽤的⽅法，但是它不适⽤于⾼维空间。

重要采样

与拒绝采样的情形相同，重要采样基于的是对提议分布q(z)的使⽤，我们很容易从提议分布中采样。之后，我们可以通过q(z)中的样本{$z^{(l)}$}的有限和的形式来表⽰期望：
$$
\begin{align}
E(f)&=\int f(z)p(z)dz
\\
&=\int f(z)\frac{p(z)}{q(z)}q(z)dz
\\
&\approx\frac{1}{L}\sum^L f(z^{(l)})\frac{p(^{(l)})}{q(^{(l)})}
\end{align}
$$
其中$r_l=\frac{p(^{(l)})}{q(^{(l)})}$​被称为重要性权重（importance weights），修正了由于从错误的概率分布中采样引⼊的偏差。注意，与拒绝采样不同，所有⽣成的样本都被保留。

常见的情形是，概率分布p(z)的计算结果没有归⼀化，即$p(z)=\tilde p(z)/Z_p$，但$Z_p$未知。类似地，我们可能希望使⽤重要采样分布$q(z)=\tilde q(z)/Z_q$，它具有相同的性质。于是我们有：
$$
\begin{align}
E(f)&=\int f(z)p(z)dz
\\
&=\frac{Z_q}{Z_p}\int f(z)\frac{\tilde p(z)}{\tilde q(z)}q(z)dz
\\
&\approx\frac{1}{L}\sum^L f(z^{(l)})\tilde r_l
\end{align}
$$
我们可以计算比值Zp/Zq：
$$
\frac{Z_p}{Z_q}=\frac{1}{Z_q}\int \tilde p(z)dz=\int \frac{\tilde p(z)}{\tilde q(z)}q(z)dz\approx \frac{1}{L}\sum \tilde r_l
$$
故
$$
E(f)\approx\sum w_lf(z^{(l)})
$$
其中
$$
w_l=\frac{\tilde r_l}{\sum_m \tilde r_m}=\frac{\frac{\tilde p^{(l)}}{q(z^{(l)})}}{\frac{\sum_m \tilde p(z^{(m)})}{q(z^{(m)})}}
$$

伪代码：

Function Importance_Sampling_Normalized(N, f, tilde_p, q_sampler, tilde_q):
  // 初始化
  sum_weighted_f_numerator = 0.0
  sum_weights_denominator = 0.0

  // 抽样和计算权重阶段
  For i from 1 to N:
    1. // 从提议分布 q(x) 中抽取一个样本
       x_i = q_sampler()

    2. // 计算该样本在目标分布下的未归一化值
       tilde_p_value_at_x_i = tilde_p(x_i)

    3. // 计算该样本在提议分布下的未归一化值
       tilde_q_value_at_x_i = tilde_q(x_i)

    4. // 安全检查：避免除以零
       If tilde_q_value_at_x_i == 0:
         // 如果 tilde_p_value_at_x_i 也为零，则此样本贡献为0，可以跳过或赋权重0
         // 如果 tilde_p_value_at_x_i 不为零，这是一个严重问题，说明 q(x) 在 p(x) 有值的地方取了0
         // 实践中需要更稳健的错误处理或确保 q 的选择是合适的
         raw_weight_i = 0.0 // 简单处理：赋权重0
         // 也可以打印警告信息: print "Warning: tilde_q(x_i) is zero!"
       Else:
         // 计算原始（未归一化的）重要性权重
         raw_weight_i = tilde_p_value_at_x_i / tilde_q_value_at_x_i
       End If

    5. // 累加分子项：f(x_i) * 原始权重
       sum_weighted_f_numerator = sum_weighted_f_numerator + f(x_i) * raw_weight_i

    6. // 累加分母项：原始权重之和
       sum_weights_denominator = sum_weights_denominator + raw_weight_i
  End For

  // 计算最终估计值
  If sum_weights_denominator == 0:
    // 如果所有权重都为零，可能表示 p 和 q 的重叠非常差，或者 N 太小
    // 或者 f(x_i) * raw_weight_i 始终为0
    print "Warning: Sum of weights is zero. Estimate might be unreliable."
    estimate_E_p_f = 0.0 // 或返回 NaN 或根据具体情况处理
  Else:
    estimate_E_p_f = sum_weighted_f_numerator / sum_weights_denominator
  End If

  Return estimate_E_p_f
End Function

MCMC

基本方法

与拒绝采样和重要采样相同，我们再⼀次从提议分布中采样。但是这次我们记录下当前状态$z^{(\tau)}$，以及依赖于这个当前状态的提议分布$z^{(\tau)}$，从⽽样本序列$z^{(1)}$,$z^{(2)}$, . . .组成了⼀个马尔科夫链。与之前⼀样，如果我们有$p(z)=\frac{\tilde p(z)}{Z_p}$，那么我们会假定对于任意的值z都可以计算$\tilde p(z)$，虽然$Z_p$的值可能位置。提议分布本⾝被选择为⾜够简单，从⽽直接采样很容易。在算法的每次迭代中，我们从提议分布中⽣成⼀个候选样本$z^{*}$，然后根据⼀个恰当的准则接受这个样本。

在基本的Metropolis算法中，我们假定提议分布是对称的，即$q(z_A|z_B)=q(z_B|z_A)$，则候选的样本被接受的概率为：
$$
A(z^{\ast},z^{(\tau)})=min(1,\frac{\tilde p(z^{\ast})}{\tilde p(z^{(\tau)})})
$$
若候选样本被接受，则$z^{(\tau+1)}=z^{*}$，否则候选样本点被丢弃，再设置$z^{\tau+1}=z^{\tau}$，然后从概率分布$q(z|z^{(\tau+1)})$​中再次抽取⼀个候选样本。与拒绝采样不同，那⾥拒绝的样本被简单地丢弃。在Metropolis算法中，当⼀个候选点被拒绝时，前⼀个样本点会被包含到最终的样本的列表中，从⽽产⽣了样本点的多个副本。

使⽤Metropolis算法从⼀个⾼斯分布中采样的简单例⼦，我们可以看出采样过程是连着的，这也就是链。

比如我们从二维高斯分布中采样N = multivariate_normal(mean=mu, cov=S)，代码如下：

X = np.mgrid[-3:3:0.05, -3:3:0.05]
Z = np.apply_along_axis(N.pdf, 0, X)
obs = N.rvs(500)
c = 0.1
z = obs.mean(axis=0)
samples = []
nrounds = 0
nsamples = 310
while len(samples) < nsamples:
    nrounds += 1
    q = multivariate_normal(z, cov=c)# 定义提议分布
    z_star = q.rvs()

    A = min(1, N.pdf(z_star) / N.pdf(z))
    if A > rand():
        samples.append(z_star)
        z = z_star
samples = np.c_[samples]

但这个方法会受到随机游⾛⾏为的影响。

比如我们考虑⼀个由整数组成的状态空间z，概率为：
$$
\begin{align}
p(z^{(\tau+1)}=z^{(\tau)})=0.5
\\
p(z^{(\tau+1)}=z^{(\tau)}+1)=0.5
\\
p(z^{(\tau+1)}=z^{(\tau)}-1)=0.5
\end{align}
$$
则根据对称性,$E[z^{(\tau)}]=0$。且

$$
\begin{align}
E_{\tau} [l(z^{(\tau)})^2] &= 0.5 \cdot E_{\tau-1}[l(z^{(\tau-1)})^2] + 0.25 \cdot E_{\tau-1}[l(z^{(\tau-1)}+1)^2] + 0.25 \cdot E_{\tau-1}[l(z^{(\tau-1)}-1)^2]
\\&= E_{\tau-1}[l(z^{(\tau-1)})^2] + 0.5
\end{align}
$$

故$E[(z^{(\tau)})^2]=\frac{\tau}{2}$。

可以看出随机游⾛在探索状态空间时是很低效的。所以设计马尔科夫链蒙特卡罗⽅法的⼀个中⼼⽬标就是避免随机游⾛⾏为。

马尔科夫链

我们可以按照下面的方式具体化一个马尔可夫链：给定初始变量的概率分布 $p(z^{(0)})$，以及后续变量的条件概率，用转移概率 (transition probability) $T_m(z^{(m)}, z^{(m+1)}) \equiv p(z^{(m+1)} | z^{(m)})$ 的形式表示。

如果对于所有的 $m$，转移概率都相同，那么这个马尔可夫链被称为同质的 (homogeneous)。

对于一个特定的变量，边缘概率可以根据前一个变量的边缘概率用链式乘积的方式表示出来，形式为

$$p(z^{(m+1)}) = \sum _ {z^{(m)}} p(z^{(m+1)} | z^{(m)}) p(z^{(m)}) $$

对于一个概率分布来说，如果马尔可夫链中的每一步都让这个概率分布保持不变，那么我们说这个概率分布关于这个马尔可夫链是不变的，或者静止的。因此，对于一个转移概率为 $T(z’, z)$ 的同质的马尔可夫链来说，如果

$$\sum _ {z’} T(z’, z) p^\ast(z’) = p^\ast(z) $$

那么概率分布 $p^\ast(z)$ 是不变的。

确保所求的概率分布 $p(z)$ 不变的一个充分（非必要）条件是令转移概率满足细节平衡 (detailed balance) 性质，定义为

$$p^\ast(z) T(z, z’) = p^\ast(z’) T(z’, z)$$

对特定的概率分布 $p^\ast(z)$ 成立。

很容易看到，满足关于特定概率分布的细节平衡性质的转移概率会使得那个概率分布具有不变性，因为

$$\sum _ {z’} p^\ast(z’) T(z’, z) = \sum _ {z’} p^\ast(z) T(z, z’) = p^\ast(z) \sum _ {z’} T(z, z’) = p^\ast(z)$$

满足细节平衡性质的马尔可夫链被称为可翻转的 (reversible)。

我们的目标是使用马尔可夫链从一个给定的概率分布中采样。如果我们构造一个马尔可夫链使得所求的概率分布是不变的，那么我们就可以达到这个目标。然而，我们还要求对于 $m \to \infty$，概率分布 $p(z^{(m)})$ 收敛于所求的不变的概率分布 $p^\ast(z)$，与初始概率分布 $p(z^{(0)})$ 无关，这种性质被称为各态历经性 (ergodicity)。这个不变的概率分布被称为均衡 (equilibrium) 分布。很明显，一个具有各态历经性的马尔可夫链只能有唯一的一个均衡分布。可以证明，同质的马尔可夫链具有各态历经性，只需对不变的概率分布和转移概率做出较弱的限制即可。

Metropolis-Hastings算法

这种情形下，提议分布不再是参数的⼀个对称函数。
$$
A_k(z^{\ast},z^{(\tau)})=min(1,\frac{\tilde p(z^{\ast})q_k(z^{(\tau)}|z^{\ast})}{\tilde p(z^{(\tau)})q_k(z^{\ast}|z^{\tau})})
$$

我们可以证明该定义的马尔科夫链是⼀个不变的概率分布，我们可以先证细节平衡，再利用充分性：
$$
\begin{align}
p(z)q_k(z’|z)A_k(z’,z)&=min(p(z)q_k(z’|z),p(z’)q_k(z|z’))
\\
&=min(p(z’)q_k(z|z’),p(z)q_k(z’|z))
\\
&=p(z’)q_k(z|z’)A_k(z,z’)
\end{align}
$$
然而，如果概率分布在不同的⽅向上的差异⾮常⼤，那么Metropolis-Hastings算法的收敛速度会⾮常慢。故有了以下改进。

吉布斯采样

例如，对于一个三维随机向量 $(x_1, x_2, x_3)$，吉布斯采样的迭代步骤如下：

1. 从 $p(x_1 | x_2^{(t)}, x_3^{(t)})$ 中抽取 $x_1^{(t+1)}$
2. 从 $p(x_2 | x_1^{(t+1)}, x_3^{(t)})$ 中抽取 $x_2^{(t+1)}$
3. 从 $p(x_3 | x_1^{(t+1)}, x_2^{(t+1)})$ 中抽取 $x_3^{(t+1)}$

为了证明这个步骤能够从所需的概率分布中采样，我们⾸先注意到对于吉布斯采样的每个步骤来说，概率分布p(z)是不变的，因此对于整个马尔科夫链来说也是不变的。这是由于当我们从$p(z_i|z _ {\backslash i})$中采样时，边缘概率分布$p(z _ {\backslash i})$显然是不变的，因为$z _ {\backslash i }$的值是不变的。并且，根据定义，对于每个步骤中来⾃正确条件概率分布$p(z_i|z _ {\backslash i})$的样本，条件概率分布都是不变的。由于条件概率分布和边缘概率分布共同确定的联合概率分布，因此我们看到联合概率分布本⾝是不变的。

为了让吉布斯采样能够从正确的概率分布中得到样本，第⼆个需要满⾜的要求为各态历经性。各态历经性的⼀个充分条件是没有条件概率分布处处为零。如果这个要求满⾜，那么z空间中的任意⼀点都可以从其他的任意⼀点经过有限步骤达到，这些步骤中每次对⼀个变量进⾏更新。如果这个要求没有满⾜，即某些条件概率分布为零，那么在这种情况下应⽤吉布斯采样时，必须显式地证明各态历经性。

吉布斯采样是MH算法的一个特例。考虑一个Metropolis-Hastings采样的步骤，它涉及到变量 $z_k$，同时保持剩余的变量 $z _ {\setminus k}$ 不变，并且对于这种情形来说，从 $z$ 到 $z^\ast$ 的转移概率为 $q _ k(z^\ast | z) = p(z^\ast_k | z _ {\setminus k})$。我们注意到 $z^\ast _ {\setminus k} = z _ {\setminus k}$，因为在采样的步骤中，向量的各个元素都不改变。并且，$p(z) = p(z_k | z _ {\setminus k})p(z _ {\setminus k})$。因此，确定Metropolis-Hastings算法中的接受概率的因子为

$$
A(z^\ast, z) = \frac{p(z^\ast)q_k(z | z^\ast)}{p(z)q_k(z^\ast | z)} = \frac{p(z^\ast_k | z^\ast _ {\setminus k})p(z^\ast _ {\setminus k})p(z_k | z^\ast _ {\setminus k})}{p(z_k | z _ {\setminus k})p(z _ {\setminus k})p(z^\ast_k | z _ {\setminus k})} = 1
$$

注意我们用到了$z _ {\backslash k}^\ast=z _ {\backslash k}$。

我们也可以证明吉布斯采样是细节平衡的。
$$
\begin{align}
p^\ast(z) T(z, z’) &=p(z_1^\tau,z_2^\tau,…,z_M^\tau)p(z_j^{\tau+1}|z_{\backslash j}^\tau)
\\&=p(z_j^\tau|z_{\backslash j}^\tau)p(z_{\backslash j}^\tau)p(z_j^{\tau+1}|z_{\backslash j}^\tau)
\\&=p(z_j^\tau|z_{\backslash j}^{\tau+1})p(z_{\backslash j}^{\tau+1})p(z_j^{\tau+1}|z_{\backslash j}^{\tau+1})
\\&=p(z_j^{\tau}|z_{\backslash j}^{\tau+1})p(z_1^{\tau+1},z_2^{\tau+1},…,z_M^{\tau+1})
\\&= p^\ast(z’) T(z’, z)
\end{align}
$$
代码：

比如我们对$p(x,y)=exp(-x^2y^2)$，则：
$$
\begin{align}
X|Y=y\sim N(0,\frac{1}{2y^2})
\\
Y|X=x\sim N(0,\frac{1}{2x^2})
\end{align}
$$

x, y = 0.01, 0.01
samples = [(x, y)]

seed(31)
for _ in range(1_000):
    sigma_x = np.sqrt(1 / (2 * y ** 2))
    x_giv_y = norm(loc=0, scale=sigma_x)
    x = x_giv_y.rvs()
    samples.append([x, y])
    
    sigma_y = np.sqrt(1 / (2 * x ** 2))
    y_giv_x = norm(loc=0, scale=sigma_y)
    y = y_giv_x.rvs()
    samples.append([x, y])
    
samples = np.r_[samples]

更多

此外，还有切⽚采样、源于模拟哈密顿动⼒学的混合蒙特卡罗⽅法等，我们在这里不再赘述。

我们也可以结合神经网络来改进MCMC。比如可以使用神经网络（例如，归一化流 Normalizing Flows）来学习复杂的、自适应的提议分布，即训练一个归一化流模型来学习并逼近目标后验分布，然后利用这个训练好的流模型作为 MCMC 算法中的提议分布。像今年新出的《Can Transformers Learn Full Bayesian Inference in Context?》

总结

MCMC的核心步骤可以概括为：

1. 构造马尔可夫链：设计一个状态转移机制（通常由一个转移核或转移概率矩阵 $P$ 定义），使得该马尔可夫链的平稳分布是目标分布 $\pi(x)$。一个关键的条件是细致平稳条件（Detailed Balance Condition），即对于任意两个状态 $i$ 和 $j$，满足 $\pi(i)P(i,j) = \pi(j)P(j,i)$。满足细致平稳条件的马尔可夫链，其平稳分布就是 $\pi(x)$。
2. 从初始状态开始迭代采样：选择一个初始状态 $x^{(0)}$，然后根据构造的马尔可夫链的转移规则，不断从当前状态转移到下一个状态，生成一系列样本 $x^{(1)}, x^{(2)}, \ldots, x^{(N)}$。
3. 达到平稳分布后收集样本：由于马尔可夫链需要一段时间才能达到平稳分布，通常会舍弃初始阶段产生的一批样本（称为“老化期”或“burn-in period”）。之后产生的样本则被认为是来自目标分布 $\pi(x)$ 的近似独立同分布样本（尽管它们在链中是相关的）。
4. 利用样本进行估计：使用收集到的有效样本来进行蒙特卡罗积分，例如估计某个函数 $f(x)$ 在目标分布 $\pi(x)$ 下的期望值：$E[f(x)] \approx \frac{1}{M} \sum _ {i=1}^{M} f(x^{(i)})$，其中 $M$ 是有效样本的数量。

有多种具体的MCMC算法：

• Metropolis-Hastings (MH) 算法：这是一个非常通用和基础的MCMC算法。它通过一个“提议分布”（proposal distribution）$q(x’|x)$ 来生成候选样本 $x’$，然后根据一个接受概率 $A(x, x’)$ 来决定是否接受这个新样本。接受概率通常定义为：
  $$A(x, x’) = \min\left(1, \frac{\pi(x’)q(x|x’)}{\pi(x)q(x’|x)}\right)$$
  如果提议分布是对称的，即 $q(x’|x) = q(x|x’)$，则MH算法简化为Metropolis算法。MH算法的优点在于它只需要计算目标分布 $\pi(x)$ 的比例，而不需要知道其归一化常数。

• 吉布斯采样 (Gibbs Sampling)：这是MH算法的一个特例，特别适用于处理高维联合分布。当直接从联合分布抽样困难，但从每个变量的全条件概率分布 (full conditional probability distribution) 抽样相对容易时，吉布斯采样非常有效。其基本思想是，轮流对每个变量进行抽样，抽样时固定其他所有变量的当前值。例如，对于一个三维随机向量 $(x_1, x_2, x_3)$，吉布斯采样的迭代步骤如下：

  1. 从 $p(x_1 | x_2^{(t)}, x_3^{(t)})$ 中抽取 $x_1^{(t+1)}$
  2. 从 $p(x_2 | x_1^{(t+1)}, x_3^{(t)})$ 中抽取 $x_2^{(t+1)}$
  3. 从 $p(x_3 | x_1^{(t+1)}, x_2^{(t+1)})$ 中抽取 $x_3^{(t+1)}$
     吉布斯采样的接受概率恒为1，因此所有提议的样本都会被接受。

MCMC的收敛性诊断

评估MCMC算法是否收敛到目标平稳分布至关重要，因为基于未收敛的链进行的推断可能是错误的。常用的收敛诊断方法包括：

• 迹图 (Trace Plots)：绘制参数样本值随迭代次数变化的图形。如果链已收敛，迹图应该看起来像一个稳定的水平带状区域，没有明显的趋势或周期性模式。
• 密度图 (Density Plots)：绘制样本的经验概率密度图。对于从不同初始点开始的多条链，如果它们都已收敛，其密度图应该相似。
• 自相关图 (Autocorrelation Plots)：显示样本序列中不同滞后期的相关性。理想情况下，自相关性应随着滞后的增加而迅速下降。高自相关意味着链的混合较差，需要更多样本。
• Gelman-Rubin 统计量 (R-hat)：通过比较多条从不同初始点并行运行的马尔可夫链的链内方差和链间方差来评估收敛性。R-hat值接近1表明链已收敛。通常认为小于1.1或1.2的值表示可接受的收敛。
• 有效样本量 (Effective Sample Size, ESS)：考虑到样本间的自相关性，ESS估计了等效的独立样本数量。ESS较低表明需要运行更长的链或改进采样器。
• 接受率 (Acceptance Rate) (主要针对MH算法)：指提议样本被接受的比例。过高或过低的接受率都可能表明采样效率低下。对于某些类型的MH算法，存在理论上的最优接受率范围（例如，对于高斯提议分布和高斯目标分布，最优接受率约为0.234）。

实践中，通常会结合使用多种诊断方法来综合判断MCMC算法的收敛性。在确认收敛后，还需要丢弃“老化期”的样本，并可能对剩余样本进行“减薄”处理，以获得用于后续分析的近似独立的样本。

总而言之，MCMC是一套功能强大的统计计算工具，它为从复杂概率分布中抽样提供了可行的解决方案，极大地推动了贝叶斯统计和相关领域的发展与应用。理解其原理、掌握常用算法并关注其收敛性是有效运用MCMC方法的关键。

参考资料

[1]A Short History of Markov Chain Monte Carlo: Subjective Recollections from Incomplete Data

[2]PRML