第 12 章 MCMC 方法基础
第 10 章介绍的 Monte Carlo 方法有一个前提：我们能够从目标分布中独立抽样。但在很多统计模型中，目标分布只知道一个未归一化的形式，或者维度较高、相关性较强，无法直接抽样。MCMC（Markov chain Monte Carlo，马尔科夫链 Monte Carlo）正是为这类问题设计的。
在贝叶斯统计中，我们经常希望从后验分布
\[ p(\theta\mid y)\propto p(y\mid\theta)p(\theta) \]
中抽样。右边的似然和先验通常容易计算，但归一化常数
\[ p(y)=\int p(y\mid\theta)p(\theta)\,d\theta \]
可能很难计算。MCMC 的关键优势是：许多算法只需要后验密度的比例，而不需要知道这个归一化常数。
12.1 马尔科夫链的基本概念
马尔科夫链是一列随机变量
\[ X_0,X_1,X_2,\ldots \]
满足如下性质：给定当前状态 \(X_t\) 后，下一状态 \(X_{t+1}\) 的条件分布不再依赖更早历史，即
\[ P(X_{t+1}\in A\mid X_t,X_{t-1},\ldots,X_0) = P(X_{t+1}\in A\mid X_t). \]
MCMC 的目标是构造一条马尔科夫链，使其长期分布等于目标分布 \(\pi(x)\)。如果链已经运行足够久，后面的样本就可以近似看作来自 \(\pi(x)\)，从而用样本平均近似期望。
与独立 Monte Carlo 不同，MCMC 样本之间通常相关。因此，样本数相同并不意味着信息量相同，这也是后面要进行诊断的原因。
12.2 Metropolis–Hastings 算法
Metropolis–Hastings 算法是最基本的 MCMC 算法之一。设目标密度为 \(\pi(x)\)，可以只知道它的比例。从当前状态 \(x_t\) 出发，先从提议分布 \(q(x^\ast\mid x_t)\) 中生成候选值 \(x^\ast\)，再以一定概率接受它。
接受概率为
\[ \alpha = \min\left\{ 1, \frac{\pi(x^\ast)q(x_t\mid x^\ast)} {\pi(x_t)q(x^\ast\mid x_t)} \right\}. \]
若接受，则 \(x_{t+1}=x^\ast\)；若拒绝，则 \(x_{t+1}=x_t\)。
Metropolis–Hastings 算法
- 给定初始值 \(x_0\)。
- 对 \(t=0,1,\ldots,T-1\)：从 \(q(x^\ast\mid x_t)\) 中生成候选值 \(x^\ast\)；计算接受概率
\[ \alpha = \min\left\{ 1, \frac{\pi(x^\ast)q(x_t\mid x^\ast)} {\pi(x_t)q(x^\ast\mid x_t)} \right\}; \]
生成 \(U\sim U(0,1)\)，若 \(U\leq\alpha\)，令 \(x_{t+1}=x^\ast\)，否则令 \(x_{t+1}=x_t\)。
实际编程时，通常在对数尺度上计算接受比，以避免上溢和下溢。
12.3 随机游走 Metropolis
若提议分布为
\[ x^\ast=x_t+\varepsilon,\quad \varepsilon\sim q(\varepsilon), \]
且 \(q(\varepsilon)=q(-\varepsilon)\)，则
\[ q(x^\ast\mid x_t)=q(x_t\mid x^\ast). \]
接受概率简化为
\[ \alpha=\min\left\{1,\frac{\pi(x^\ast)}{\pi(x_t)}\right\}. \]
这就是随机游走 Metropolis 算法。提议步长很重要：步长太小，链移动慢；步长太大，候选值经常被拒绝。
下面从一个双峰混合正态分布中抽样。这个例子展示 MCMC 可以只利用未归一化密度的相对大小。
set.seed(2026)
target <- function(x) {
0.35 * dnorm(x, mean = -2, sd = 0.7) +
0.65 * dnorm(x, mean = 2, sd = 1)
}
rw_metropolis <- function(x0, n_iter, proposal_sd) {
x <- numeric(n_iter)
x[1] <- x0
accept <- 0
for (t in 2:n_iter) {
proposal <- rnorm(1, mean = x[t - 1], sd = proposal_sd)
log_alpha <- log(target(proposal)) - log(target(x[t - 1]))
if (log(runif(1)) < log_alpha) {
x[t] <- proposal
accept <- accept + 1
} else {
x[t] <- x[t - 1]
}
}
list(samples = x, accept_rate = accept / (n_iter - 1))
}
fit <- rw_metropolis(x0 = 0, n_iter = 10000, proposal_sd = 1.2)
c(accept_rate = fit$accept_rate,
sample_mean = mean(fit$samples[2001:10000]))
#> accept_rate sample_mean
#> 0.6647 0.7297这里前 2000 次迭代被视为 burn-in，即为了减弱初始值影响而舍弃的早期样本。burn-in 的长度没有统一公式，需要结合 trace plot 和具体问题判断。
12.4 Gibbs 抽样
Gibbs 抽样是另一类常见 MCMC 方法。若多维目标分布为 \(\pi(x_1,\ldots,x_d)\)，并且每个条件分布
\[ \pi(x_j\mid x_{-j}) \]
都容易抽样，就可以轮流从条件分布中更新每个分量。这里 \(x_{-j}\) 表示除 \(x_j\) 以外的其他分量。
以二维正态分布为例，若
\[ \begin{pmatrix}X\\Y\end{pmatrix} \sim N\left[ \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix}1&\rho\\ \rho&1\end{pmatrix} \right], \]
则条件分布为
\[ X\mid Y=y\sim N(\rho y,1-\rho^2), \]
\[ Y\mid X=x\sim N(\rho x,1-\rho^2). \]
set.seed(1)
rho <- 0.8
n_iter <- 5000
xy <- matrix(0, n_iter, 2)
for (t in 2:n_iter) {
xy[t, 1] <- rnorm(1, mean = rho * xy[t - 1, 2],
sd = sqrt(1 - rho^2))
xy[t, 2] <- rnorm(1, mean = rho * xy[t, 1],
sd = sqrt(1 - rho^2))
}
colMeans(xy[1001:n_iter, ])
#> [1] -0.019701 -0.008747
cor(xy[1001:n_iter, ])
#> [,1] [,2]
#> [1,] 1.0000 0.8059
#> [2,] 0.8059 1.0000Gibbs 抽样在层次模型、缺失数据模型和一些贝叶斯回归模型中非常有用。它的限制是：条件分布必须容易抽样；若变量之间相关性很强，链也可能移动缓慢。
12.5 案例：信用违约概率的贝叶斯逻辑回归
下面构造一个模拟的信用风险数据。令 \(y_i=1\) 表示违约，\(x_i\) 表示标准化收入。逻辑回归模型为
\[ P(y_i=1\mid x_i,\beta) = \frac{\exp(\beta_0+\beta_1x_i)} {1+\exp(\beta_0+\beta_1x_i)}. \]
给定正态先验
\[ \beta_0,\beta_1\sim N(0,10^2), \]
后验分布没有简单闭式形式。我们用随机游走 Metropolis 从\(p(\beta_0,\beta_1\mid y,x)\) 中抽样。数据仍然是模拟数据。
set.seed(2)
n <- 300
income <- as.numeric(scale(rlnorm(n, meanlog = log(8), sdlog = 0.55)))
eta <- -1.0 - 1.1 * income
default <- rbinom(n, size = 1, prob = plogis(eta))
log_posterior <- function(beta) {
eta <- beta[1] + beta[2] * income
loglik <- sum(default * eta - log1p(exp(eta)))
logprior <- sum(dnorm(beta, mean = 0, sd = 10, log = TRUE))
loglik + logprior
}
n_iter <- 8000
samples <- matrix(NA, n_iter, 2)
samples[1, ] <- c(0, 0)
accept <- 0
proposal_sd <- c(0.12, 0.12)
for (t in 2:n_iter) {
current <- samples[t - 1, ]
proposal <- current + rnorm(2, sd = proposal_sd)
log_alpha <- log_posterior(proposal) - log_posterior(current)
if (log(runif(1)) < log_alpha) {
samples[t, ] <- proposal
accept <- accept + 1
} else {
samples[t, ] <- current
}
}
post <- samples[2001:n_iter, ]
colnames(post) <- c("beta0", "beta1")
rbind(mean = colMeans(post),
q025 = apply(post, 2, quantile, 0.025),
q975 = apply(post, 2, quantile, 0.975))
#> beta0 beta1
#> mean -0.8648 -1.0737
#> q025 -1.1640 -1.4722
#> q975 -0.5893 -0.7083
accept / (n_iter - 1)
#> [1] 0.6281后验均值和区间反映了参数不确定性。若 \(\beta_1\) 的后验大部分小于 0，可以解释为：在这个模拟模型中，收入越高，违约概率越低。真实信贷分析还需要更多变量、更严格的数据来源说明和模型诊断。
12.6 MCMC 诊断
MCMC 输出不能只看后验均值。至少应检查以下方面：
- trace plot：样本轨迹是否在合理范围内来回移动，而不是卡在某个区域。
- 接受率：随机游走 Metropolis 接受率过高可能表示步长太小，过低可能表示步长太大。
- 自相关：相邻样本相关性越强，有效信息量越少。
- 多链比较：从不同初始值出发，检查是否得到相似后验结果。
- 有效样本量：关注独立信息量，而不是只看迭代次数。
这些诊断不是形式要求，而是判断模拟近似是否可信的必要步骤。第 20 章会介绍现代贝叶斯软件中常见的 R-hat 和有效样本量指标。