PricingLoginStart free
R Academy · 课时

MCMC 抽样与诊断

运行抽样、检查链，并解读 Rhat 和 ESS 诊断结果

第 3 / 4 课13 个步骤

MCMC 抽样与诊断 是 CoddyKit 上的免费 R Academy 课时。 这是第 3 节课，共 4 节。 你可以在下方免费阅读本课时的完整内容 — 然后在浏览器中使用内置代码编辑器和全天候 AI 导师进行实践。 这是 R Academy 学习路径的一部分，你的进度在网页和 CoddyKit 应用中同步。 R Academy 课程共包含 4 节课。

什么是 MCMC

马尔可夫链蒙特卡罗（MCMC）是一类算法，用于在无法直接抽样时从概率分布中抽取样本。在贝叶斯统计中，MCMC 从后验分布 P(parameters | data) 中抽样。

RStan 实现了无 U 形回转抽样器（NUTS），这是一种先进的 MCMC 算法。

一个简单的 Stan 模型

Stan 模型是一个文本块，用于定义数据类型、参数和对数后验。最简单的模型是在方差已知的情况下估计正态分布的均值。

# library(rstan)
#
# stan_code <- '
# data {
#   int<lower=0> N;
#   vector[N] y;
# }
# parameters {
#   real mu;
#   real<lower=0> sigma;
# }
# model {
#   mu    ~ normal(0, 10);   // prior
#   sigma ~ exponential(1);   // prior
#   y     ~ normal(mu, sigma); // likelihood
# }
# '

调用 stan() 进行抽样

stan(model_code=..., data=..., chains=4, iter=2000, warmup=1000) 会先编译模型（仅一次），然后抽取样本。使用 iter=2000 和 warmup=1000 时，每条链会产生 1000 个预热后的样本；4 条链总共产生 4000 个样本。

# library(rstan)
# options(mc.cores = parallel::detectCores())
#
# y <- c(2.1, 1.8, 2.4, 1.9, 2.3, 2.0, 1.7, 2.2)
# stan_data <- list(N = length(y), y = y)
#
# fit <- stan(
#   model_code = stan_code,
#   data       = stan_data,
#   chains     = 4,
#   iter       = 2000,
#   warmup     = 1000,
#   seed       = 42
# )

print(fit) — 汇总表

print(fit) 会为每个参数显示一张汇总表，其中包括后验均值、标准差、分位数、Rhat 和 n_eff。这两个诊断指标是首先需要检查的内容。

# print(fit)
#
# Example output:
#       mean se_mean   sd  2.5%   25%   50%   75%  97.5%  n_eff Rhat
# mu    2.05    0.00 0.15  1.76  1.95  2.05  2.15   2.34   3842    1
# sigma 0.22    0.00 0.06  0.13  0.18  0.21  0.25   0.37   3521    1
# lp__  4.38    0.02 1.01  1.60  3.91  4.71  5.18   5.50   2148    1

Rhat 收敛标准

Rhat（潜在尺度缩减因子）会比较链内方差与链间方差。接近 1.0 的值表明所有链都已收敛到同一个分布。

  • Rhat < 1.01 — 已收敛（当前标准）
  • Rhat > 1.01 — 链尚未充分混合；请运行更多迭代
  • Rhat > 1.1 — 存在严重的收敛问题
# Check Rhat for all parameters:
# s <- summary(fit)$summary
# rhat_vals <- s[, 'Rhat']
# cat('Max Rhat:', max(rhat_vals, na.rm = TRUE), '
')
# if (any(rhat_vals > 1.01, na.rm = TRUE)) {
#   warning('Convergence issue detected!')
# } else {
#   cat('All Rhat < 1.01 — chains converged
')
# }

n_eff — 有效样本量

n_eff（有效样本量）会考虑连续 MCMC 样本之间的自相关。相关样本包含的信息少于独立样本。

  • n_eff 接近总迭代次数 — 样本近似独立，效果 excellent
  • n_eff / total_samples > 0.1 — 通常可以接受
  • n_eff 很低 — 自相关很高；请考虑重新参数化模型
# s <- summary(fit)$summary
# n_eff_vals <- s[, 'n_eff']
# total_samples <- 4 * 1000   # chains * post-warmup iter
# ratio <- n_eff_vals / total_samples
# cat('n_eff ratio (mu)   :', round(ratio['mu'], 2), '
')
# cat('n_eff ratio (sigma):', round(ratio['sigma'], 2), '
')

使用 traceplot() 进行可视化收敛检查

traceplot(fit, pars = 'mu') 会绘制每条链在各次迭代中抽取的 mu 值。已经收敛的链看起来像一条模糊的毛毛虫：所有链彼此重叠，没有趋势或漂移。未收敛的链会四处游走或保持分离。

# library(rstan)
#
# traceplot(fit, pars = c('mu', 'sigma'), inc_warmup = FALSE)
#
# Good traceplot characteristics:
# - All 4 chains overlapping completely (same range)
# - No visible drift or trend
# - Rapid mixing (values jump around quickly)
# - No flat regions (stuck sampler)
cat('A healthy traceplot looks like a fuzzy caterpillar
')

使用 pairs() 查看后验相关性

pairs(fit, pars = c('mu', 'sigma')) 会显示后验样本的散点图矩阵。它可以揭示参数之间的相关性，并突出显示发散跃迁（以红色绘制）；这些跃迁表明抽样器难以处理相应区域。

# pairs(fit, pars = c('mu', 'sigma'))
#
# What to look for:
# - Elliptical clouds: mild correlation (OK)
# - Banana / funnel shapes: reparameterization needed
# - Red dots (divergences): geometry problem in posterior
#   => increase adapt_delta: stan(..., control=list(adapt_delta=0.95))
cat('Red dots in pairs() indicate divergent transitions — investigate!
')

提取后验样本

extract(fit, pars = 'mu')$mu 会返回 mu 的所有预热后样本组成的数值向量。您可以使用这些样本计算任意后验汇总，例如均值、可信区间和某个条件成立的概率。

# mu_samples <- extract(fit, pars = 'mu')$mu
# cat('Posterior mean :', mean(mu_samples), '
')
# cat('95% CI:', quantile(mu_samples, c(0.025, 0.975)), '
')
# cat('P(mu > 2):', mean(mu_samples > 2), '
')
# hist(mu_samples, main = 'Posterior of mu', xlab = 'mu', col = 'steelblue')

启动 ShinyStan 进行交互式诊断

shinystan::launch_shinystan(fit) 会打开一个交互式 Shiny 应用，在同一处提供轨迹图、后验分布图、成对图和 NUTS 诊断。它是探索 RStan 拟合结果最全面的工具。

# install.packages('shinystan')
# library(shinystan)
#
# shinystan::launch_shinystan(fit)
#
# ShinyStan tabs:
# - Diagnose: Rhat, n_eff, divergences, energy
# - Explore:  marginal posteriors, scatter plots
# - Model:    Stan code, data
# - NUTS:     step size, tree depth per chain

常见的收敛问题修复方法

当 Rhat > 1.01 或您看到发散时：

  • 增加 iter 和 warmup
  • 在 control 中将 adapt_delta 提高到接近 1.0（例如 0.95）
  • 重新参数化 — 对层级模型使用非中心化参数化
  • 如果先验过于分散，请收紧先验
  • 检查数据 — 异常值或尺度差异会导致抽样器出现问题
# Re-run with higher adapt_delta to reduce divergences:
# fit2 <- stan(
#   model_code = stan_code,
#   data       = stan_data,
#   chains     = 4,
#   iter       = 4000,
#   warmup     = 2000,
#   control    = list(adapt_delta = 0.95, max_treedepth = 12),
#   seed       = 42
# )

快速检查：Rhat 阈值

表示 Stan 模型已经收敛的 Rhat 当前标准阈值是多少？

在交互式课程中回答

MCMC 抽样与诊断回顾

RStan MCMC 的主要工作流程：

  • stan(model_code=..., data=..., chains=4, iter=2000, warmup=1000) 拟合模型
  • print(fit) 显示 Rhat 和 n_eff — 主要的收敛诊断指标
  • Rhat < 1.01 且 n_eff / total > 0.1 表明样本表现良好
  • traceplot() — 可视化混合检查；pairs() — 揭示后验几何结构问题
  • extract(fit, pars='mu')$mu — 访问原始后验样本
  • shinystan::launch_shinystan(fit) — 全面的交互式诊断
免费开始

用 AI 导师学习 R — 免费

在浏览器中编写并运行真实代码，获得全天候 AI 导师的即时帮助，并在网页或应用中继续学习。

课程
43
课程
159
免费开始学习查看完整课程

常见问题解答

「MCMC 抽样与诊断」课时是免费的吗？

是的 — 「MCMC 抽样与诊断」的完整文本可在网页上免费阅读。要进行交互式练习（内置代码编辑器和全天候 AI 导师）并解锁 R Academy 课程的其余内容，请升级到 CoddyKit PRO。 R Academy 课程共包含 4 节课。

「MCMC 抽样与诊断」这节课中我会学到什么？

运行抽样、检查链，并解读 Rhat 和 ESS 诊断结果 你通过在浏览器中直接运行的动手代码来练习 R Academy，全天候 AI 导师会在你学习这节课的过程中回答你的问题。

学习 R Academy 需要有经验吗？

无需任何先前经验。CoddyKit 上的 R Academy 课程适合初学者到高级学习者，你可以从这里开始或从头开始，按照自己的节奏学习。 这是第 3 节课，共 4 节。

「MCMC 抽样与诊断」课时需要多长时间？

大多数 CoddyKit 课程大约需要 5–10 分钟。每节课都很精短且互动，所以你能稳步进步，并在网页和应用中从离开的地方继续。

我能在这节 R Academy 课中编写并运行代码吗？

能。每节 R Academy 课都包含内置代码编辑器，你可以在浏览器中直接编写并运行真实代码，并获得即时 AI 反馈 — 无需本地设置。

此课程中的所有课时

  1. 贝叶斯思维入门
  2. 在 R 中编写 Stan 模型
  3. MCMC 抽样与诊断
  4. 后验预测检验
← 返回 R Academy