在bugs/stan模型中提供不同的先验族作为参数

Question

这是Andrew Gelman在贝叶斯数据分析中的经典八派例子。请参阅下面的标准文件和R代码。我在stan文件中使用了一个带有对聚体a的cauchy先验作为超参数tau。我正试图提供R函数"；学校；具有不在cauchy族中的不同先验，例如，一致(01000(先验，这样我就不必为新的先验创建不同的stans文件。这可能发生在stan或bug中吗？

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);
}

`

`

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()

`

Answer 1

可以在Stan代码中预先指定多个分布，然后在输入数据中指定您想要的分布。斯坦并不是真的打算用这种方式，但它是可以做到的！

下面是一个例子。我添加了一个新的数据变量tau_prior；它是一个整数，用于指定要用于tau的哪个先验。1=柯西，2=均匀，3=指数。此外，对于每种类型的先验，都有一个数据变量来设置超参数。(未选择的分布的超参数没有影响。(

data {
int<lower=0> J;         // number of schools
real y[J];              // estimated treatment effects
real<lower=0> sigma[J]; // standard error of effect estimates
int<lower=1,upper=3> tau_prior;
real<lower=0> cauchy_sigma;
real<lower=0> uniform_beta;
real<lower=0> exponential_beta;
}
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);
if(tau_prior == 1) {
tau ~ cauchy(0, cauchy_sigma);
} else if(tau_prior == 2) {
tau ~ uniform(0, uniform_beta);
} else if(tau_prior == 3) {
tau ~ exponential(exponential_beta);
}
}

我还修改了R函数，使其为每个超参数提供默认值，其范围与您已经使用的范围相似。

school <- function(tau_prior = 1,
cauchy_sigma = 100,
uniform_beta = 1000,
exponential_beta = 0.01) {
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),
tau_prior = tau_prior,
cauchy_sigma = cauchy_sigma,
uniform_beta = uniform_beta,
exponential_beta = exponential_beta)
fit <- stan(file = "schools.stan", data = schools_dat, iter = 20)
print(fit)
}
# The default: use a Cauchy prior with scale 100.
school()
# Use a uniform prior with the default upper limit (1000).
school(tau_prior = 2)
# Use an exponential prior with a non-default rate (1).
school(tau_prior = 3, exponential_beta = 1)