PricingLoginStart free
R Academy · 课时

在 R 中编写 Stan 模型

使用 Stan 语法定义数据块、参数和模型块

第 2 / 4 课13 个步骤

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

什么是 Stan

Stan 是一种用于贝叶斯统计建模的概率编程语言。RStan 是 Stan 的 R 接口。您使用 Stan 语言（类似 C++）编写模型，Stan 会将其编译为高效的 C++ 代码，并运行哈密顿蒙特卡罗（HMC）抽样来近似后验。

# Stan installation check
library(rstan)

# Check version
cat('RStan version:', as.character(packageVersion('rstan')), '\n')

# Enable parallel chains
options(mc.cores = parallel::detectCores())

# Reuse compiled models across sessions
rstan_options(auto_write = TRUE)

cat('Stan ready. Cores:', parallel::detectCores(), '\n')

Stan 模型结构

Stan 模型最多包含六个命名代码块：data、transformed data、parameters、transformed parameters、model 和 generated quantities。其中三个必需代码块是 data、parameters 和 model。

library(rstan)

# Stan model as a character string in R
stan_code <- '
data {
  int<lower=0> N;       // number of observations
  vector[N] x;          // predictor
  vector[N] y;          // response
}
parameters {
  real alpha;           // intercept
  real beta;            // slope
  real<lower=0> sigma;  // noise (must be positive)
}
model {
  // Priors
  alpha ~ normal(0, 10);
  beta  ~ normal(0, 10);
  sigma ~ exponential(1);

  // Likelihood
  y ~ normal(alpha + beta * x, sigma);
}
'

cat('Stan model defined as a string in R\n')
cat('Blocks: data, parameters, model\n')

data{}：声明输入数据

data 代码块声明模型从 R 接收的所有外部数据。类型包括 int、real、vector[N]、matrix[M,N] 和 array。<lower=0> 等约束会在运行时进行检查。

library(rstan)

# Comprehensive data block examples
data_block_examples <- '
data {
  // Scalars
  int<lower=1> N;               // at least 1 observation
  int<lower=2> K;               // at least 2 groups

  // Constrained scalars
  real<lower=0, upper=1> rate;  // probability

  // Vectors
  vector[N] y;                  // continuous response
  array[N] int<lower=0, upper=1> z;  // binary outcomes

  // Matrix
  matrix[N, K] X;               // design matrix

  // Integer array
  array[N] int group;           // group membership 1..K
}
'

cat(data_block_examples)

parameters{}：参数类型

parameters 块声明 Stan 将进行抽样的未知量。此块中的约束定义了参数空间：<lower=0> 表示正值，<lower=0, upper=1> 表示概率。对于受约束的参数，Stan 会自动应用对数概率的雅可比修正。

library(rstan)

# Common parameter declarations
param_examples <- '
parameters {
  // Unconstrained
  real mu;                        // mean
  vector[K] beta;                 // regression coefficients

  // Positive (sigma, lambda, variance)
  real<lower=0> sigma;
  real<lower=0> lambda;

  // Probability
  real<lower=0, upper=1> theta;

  // Simplex (sums to 1, for mixture weights)
  simplex[K] pi;

  // Correlation matrix
  corr_matrix[K] Omega;

  // Cholesky factor of covariance
  cholesky_factor_cov[K] L_Sigma;
}
'
cat(param_examples)

model{}：先验与似然

model 块通过 ~（波浪号）语法累积对数后验。y ~ normal(mu, sigma) 是将正态 PDF 的对数加到 target 上的简写形式。您也可以使用 target += normal_lpdf(y | mu, sigma) 显式写出对数概率增量。

library(rstan)

model_block_example <- '
model {
  // --- Priors ---
  mu    ~ normal(0, 10);          // weakly informative
  sigma ~ cauchy(0, 2.5);         // half-Cauchy for scale
  beta  ~ normal(0, 1);           // standardised coefficients

  // --- Likelihood ---
  // Tilde notation (most common)
  y ~ normal(mu + X * beta, sigma);

  // Equivalent explicit notation:
  // target += normal_lpdf(y | mu + X * beta, sigma);

  // For loop (less common but valid)
  // for (i in 1:N)
  //   target += normal_lpdf(y[i] | mu, sigma);
}
'
cat(model_block_example)

拟合一个简单的正态模型

下面是一个完整且最精简的 Stan 工作流程：定义模型字符串，准备数据列表，调用 stan()，并使用 print() 查看结果。Stan 第一次运行时会编译模型，并将其缓存下来供后续运行使用。

library(rstan)

# Stan model: estimate mean and SD of a normal distribution
normal_model <- '
data {
  int<lower=0> N;
  vector[N] y;
}
parameters {
  real mu;
  real<lower=0> sigma;
}
model {
  mu    ~ normal(0, 10);
  sigma ~ exponential(0.1);
  y     ~ normal(mu, sigma);
}
'

# Simulate data
set.seed(42)
y_data <- rnorm(50, mean = 5, sd = 2)

# Fit the model
fit <- stan(
  model_code = normal_model,
  data       = list(N = length(y_data), y = y_data),
  chains     = 2,
  iter       = 1000,
  warmup     = 500,
  refresh    = 0  # suppress iteration output
)

print(fit, pars = c('mu', 'sigma'))

在 Stan 中进行线性回归

Stan 让贝叶斯线性回归变得十分直接：为系数指定正态先验，为 sigma 指定指数先验或半柯西先验。后验会给出每个系数完整的不确定性分布，而不仅仅是点估计。

library(rstan)

lin_reg_model <- '
data {
  int<lower=0> N;
  vector[N] x;
  vector[N] y;
}
parameters {
  real alpha;
  real beta;
  real<lower=0> sigma;
}
model {
  alpha ~ normal(0, 10);
  beta  ~ normal(0, 10);
  sigma ~ exponential(1);
  y ~ normal(alpha + beta * x, sigma);
}
generated quantities {
  vector[N] y_rep;  // posterior predictive
  for (i in 1:N)
    y_rep[i] = normal_rng(alpha + beta * x[i], sigma);
}
'

set.seed(7)
n <- 80
x <- rnorm(n); y <- 2 + 3 * x + rnorm(n, 0, 1)

fit <- stan(model_code = lin_reg_model,
            data = list(N = n, x = x, y = y),
            chains = 2, iter = 1000, refresh = 0)
print(fit, pars = c('alpha', 'beta', 'sigma'))

transformed parameters{} 块

transformed parameters 块根据抽样得到的参数计算派生量。这些派生量有助于以数值稳定的方式参数化模型（例如使用 Cholesky 分解和对数变换），同时保持主要参数的可解释性。

library(rstan)

# Model with transformed parameters
trans_param_example <- '
data {
  int<lower=0> N;
  array[N] int<lower=0> y;  // counts
}
parameters {
  real log_lambda;           // work in log space for stability
}
transformed parameters {
  real<lower=0> lambda;
  lambda = exp(log_lambda);  // transform back to original scale
}
model {
  log_lambda ~ normal(1, 2);  // prior on log scale
  y ~ poisson(lambda);
}
'

set.seed(42)
y_counts <- rpois(30, lambda = 5)

fit_poisson <- stan(
  model_code = trans_param_example,
  data       = list(N = length(y_counts), y = y_counts),
  chains     = 2, iter = 1000, refresh = 0
)
print(fit_poisson, pars = c('log_lambda', 'lambda'))

generated quantities{} 块

generated quantities 块在抽样完成后运行，用于计算其他量：后验预测样本（y_rep）、LOO-CV 的对数似然，或用于汇总的变换参数。此处的值均从后验预测分布中抽样得到。

library(rstan)

# Use generated quantities for posterior predictive checks
model_with_gq <- '
data {
  int<lower=0> N;
  vector[N] y;
}
parameters {
  real mu;
  real<lower=0> sigma;
}
model {
  mu    ~ normal(0, 10);
  sigma ~ exponential(0.5);
  y     ~ normal(mu, sigma);
}
generated quantities {
  vector[N] y_rep;           // replicated datasets
  real mean_y_rep;           // mean of replicated data
  for (i in 1:N)
    y_rep[i] = normal_rng(mu, sigma);
  mean_y_rep = mean(y_rep);
}
'

set.seed(1)
y_obs <- rnorm(40, 3, 1.5)
fit <- stan(model_code = model_with_gq,
            data = list(N = length(y_obs), y = y_obs),
            chains = 2, iter = 1000, refresh = 0)
print(fit, pars = c('mu', 'sigma', 'mean_y_rep'))

将数据从 R 传递到 Stan

stan() 的 data 参数是一个带名称的 R 列表。名称必须与 Stan data{} 块中声明的变量名称完全匹配。向量会变成 Stan 的 vector[N]；整数会变成 int；R 矩阵会变成 Stan 的 matrix[M,N]。

library(rstan)

# Data preparation: names must match Stan data block exactly
set.seed(42)
n <- 60
x1 <- rnorm(n)
x2 <- rnorm(n)
y  <- 1.5 + 2 * x1 - 0.8 * x2 + rnorm(n, 0, 0.5)

# Build the design matrix
X <- cbind(x1, x2)  # 60 x 2 matrix

# Named list passed to stan(data = ...)
stan_data <- list(
  N = n,          # int
  K = ncol(X),    # int
  X = X,          # matrix[N, K]
  y = y           # vector[N]
)

cat('Stan data list elements:\n')
for (nm in names(stan_data)) {
  cat(' ', nm, ': class =', class(stan_data[[nm]]),
      'dim =', paste(dim(stan_data[[nm]]), collapse = 'x'),
      '\n')
}

查看已编译的模型

拟合完成后，使用 print(fit) 查看所有参数的后验汇总（mean、se_mean、sd、quantiles、Rhat、n_eff）。使用 stan_plot(fit) 绘制可视化区间图，并使用 traceplot(fit) 评估链的混合情况。

library(rstan)

# Reuse the simple normal model from scene 6
normal_model <- '
data { int<lower=0> N; vector[N] y; }
parameters { real mu; real<lower=0> sigma; }
model {
  mu ~ normal(0, 10);
  sigma ~ exponential(0.1);
  y ~ normal(mu, sigma);
}
'

set.seed(5)
fit <- stan(model_code = normal_model,
            data = list(N = 50, y = rnorm(50, 7, 3)),
            chains = 2, iter = 1000, refresh = 0)

# Detailed summary table
print(fit)

# Extract as data frame
posterior_df <- as.data.frame(fit)
cat('\nPosterior samples shape:', nrow(posterior_df),
    'rows x', ncol(posterior_df), 'cols\n')

快速检查

在 Stan 模型中，如果您希望将参数约束为严格正值（例如标准差），哪种声明是正确的？

在交互式课程中回答

回顾：编写 Stan 模型

要点：

  • Stan 模型有三个基本块：data{}、parameters{}、model{}
  • 约束：概率可使用 <lower=0>、<upper=1>、<lower=0, upper=1>
  • 类型：int、real、vector[N]、matrix[M,N]、simplex[K]
  • 模型块：使用波浪号语法 y ~ normal(mu, sigma) 表示先验和似然
  • 使用 transformed parameters{} 表示派生量；使用 generated quantities{} 表示抽样后的计算
  • 将数据作为带名称的 R 列表传入，并使其名称与 Stan 块中的名称完全匹配
  • rstan_options(auto_write = TRUE) 会缓存已编译的模型
library(rstan)

# Stan model skeleton
model_skeleton <- '
data    { int N; vector[N] y; }
parameters { real mu; real<lower=0> sigma; }
model   { mu ~ normal(0,10); sigma ~ exponential(1); y ~ normal(mu, sigma); }
'
cat(model_skeleton)
cat('\nFit with: stan(model_code = model_skeleton, data = list(N=..., y=...), chains=4)\n')
免费开始

用 AI 导师学习 R — 免费

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

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

常见问题解答

「在 R 中编写 Stan 模型」课时是免费的吗？

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

「在 R 中编写 Stan 模型」这节课中我会学到什么？

使用 Stan 语法定义数据块、参数和模型块 你通过在浏览器中直接运行的动手代码来练习 R Academy，全天候 AI 导师会在你学习这节课的过程中回答你的问题。

学习 R Academy 需要有经验吗？

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

「在 R 中编写 Stan 模型」课时需要多长时间？

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

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

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

此课程中的所有课时

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