如何在Stan/BUGS中为八校模型动态传入不同族先验？

问题 #

这是Andrew Gelman贝叶斯数据分析经典的八校示例，当前Stan文件中超参数tau使用带参数A的Cauchy先验。希望在R函数school中传入Cauchy以外的不同先验（比如uniform(0,1000)），且不需要创建多个Stan文件，请问在Stan或BUGS中是否可行？

现有schools.stan文件：

data {
  int<lower=0> J;         // number of schools 
  real y[J];              // estimated treatment effects
  real<lower=0> sigma[J]; // standard error of effect estimates 
  real<lower=0> A;
}
parameters {
  real mu;                // population treatment effect
  real<lower=0> tau;      // standard deviation in treatment effects
  vector[J] eta;          // unscaled deviation from mu by school
}
transformed parameters {
  vector[J] theta = mu + tau * eta;        // school treatment effects
}
model {
 eta ~ normal(0, 1);
y ~ normal(theta, sigma); 
tau ~ cauchy(0,A);
}

对应的R函数：

school <- function(A=100){
 schools_dat <- list(J = 8, 
                    y = c(28,  8, -3,  7, -1,  1, 18, 12),
                    sigma = c(15, 10, 16, 11,  9, 11, 10, 18),
                    A=A)

fit <- stan(file = "schools.stan", data = schools_dat,iter = 20)
print(fit)
}
school()

尝试过的代码（但不知道如何修改Stan文件）：

school <- function(prior="dunif(0,1000"){
 schools_dat <- list(J = 8, 
                    y = c(28,  8, -3,  7, -1,  1, 18, 12),
                    sigma = c(15, 10, 16, 11,  9, 11, 10, 18),
                    prior=prior)

fit <- stan(file = "schools.stan", data = schools_dat,iter = 20)
print(fit)
}
school()

解决方案 #

在Stan中实现（推荐） #

Stan不支持直接传递字符串形式的分布表达式，但可以通过数据定义先验类型标识+对应参数的方式，在单个Stan文件中支持多种先验，无需创建多个文件。具体步骤如下：

1. 修改Stan文件，增加先验类型标识变量和对应参数，在model块中用条件语句选择先验：

data {
  int<lower=0> J;         // 学校数量
  real y[J];              // 估计的处理效应
  real<lower=0> sigma[J]; // 效应估计的标准误
  int<lower=1, upper=3> prior_type; // 先验类型标识：1=Cauchy, 2=Uniform, 3=Half-Normal
  // 对应不同先验的参数
  real<lower=0> cauchy_scale; 
  real<lower=0> uniform_low;
  real<lower=uniform_low> uniform_high;
  real<lower=0> half_normal_scale;
}
parameters {
  real mu;                // 总体处理效应
  real<lower=0> tau;      // 处理效应的标准差
  vector[J] eta;          // 各学校相对于mu的未缩放偏差
}
transformed parameters {
  vector[J] theta = mu + tau * eta; // 各学校的处理效应
}
model {
  eta ~ normal(0, 1);
  y ~ normal(theta, sigma); 
  
  // 根据先验类型选择对应分布
  if (prior_type == 1) {
    tau ~ cauchy(0, cauchy_scale);
  } else if (prior_type == 2) {
    tau ~ uniform(uniform_low, uniform_high);
  } else if (prior_type == 3) {
    tau ~ normal(0, half_normal_scale) T[0, ]; // 半正态分布（截断在0以上）
  }
}

2. 修改R函数，接收先验类型和对应参数，组装成Stan所需的数据列表：

school <- function(prior_type = 1, 
                   cauchy_scale = 100, 
                   uniform_low = 0, uniform_high = 1000,
                   half_normal_scale = 50) {
  schools_dat <- list(
    J = 8, 
    y = c(28,  8, -3,  7, -1,  1, 18, 12),
    sigma = c(15, 10, 16, 11,  9, 11, 10, 18),
    prior_type = prior_type,
    cauchy_scale = cauchy_scale,
    uniform_low = uniform_low,
    uniform_high = uniform_high,
    half_normal_scale = half_normal_scale
  )
  
  fit <- stan(file = "schools.stan", data = schools_dat, iter = 2000) # 建议增加迭代次数保证收敛
  print(fit)
}

# 使用示例：
school() # 默认用Cauchy先验（scale=100）
school(prior_type = 2) # 使用Uniform(0,1000)先验
school(prior_type = 3, half_normal_scale = 30) # 使用半正态先验（scale=30）

在BUGS中实现 #

BUGS同样支持通过数据变量控制先验选择，核心思路和Stan一致：在数据中定义先验类型标识，然后在模型块中用条件语句选择对应先验分布。例如：

model {
  # 其他模型部分...
  if (prior_type == 1) {
    tau ~ dcauchy(0, cauchy_scale)
  } else if (prior_type == 2) {
    tau ~ dunif(uniform_low, uniform_high)
  }
}

对应的R调用（比如用R2OpenBUGS）时，同样传递prior_type和对应参数到数据列表即可。

注意事项 #

• 不要尝试传递字符串形式的分布表达式（比如"dunif(0,1000)"），Stan/BUGS都不会解析这类字符串作为模型代码，会直接报错。
• 增加新先验类型时，只需扩展prior_type的取值范围、补充对应数据参数和条件分支即可。
• 实际使用时建议将迭代次数从示例中的20增加到2000或更多，保证采样收敛。

内容的提问来源于stack exchange，提问作者Statisfun