Supply different families of priors as a parameter in the bugs/stan model

Question

This is the classic eight school example in Bayesian data analysis by Andrew Gelman. Please see the stan file and R code below. I use a cauchy prior with paratmer A for the hyperparamter tau in the stan file. I am trying to supply the R function "school" with different priors not within cauchy family, for example, uniform(0,1000) prior, so that I do not have to create different stans file for the new priors. Is this possible within stan or 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);
}

`

`

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

`

I tried the following but have no idea how to change the stan file correspondingly. `

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

It's possible to pre-specify more than one distribution in the Stan code, and then specify which distribution you want in the input data. Stan isn't really intended to be used this way, but it can be done!

Here's an example. I've added a new data variable, tau_prior; it's an integer that specifies which prior you want to use for tau. 1 = Cauchy, 2 = uniform, 3 = exponential. In addition, for each type of prior, there's a data variable that sets a hyperparameter. (Hyperparameters for the distributions that aren't chosen have no effect.)

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

I've also modified the R function so that it provides default values for each hyperparameter, on a scale similar to the one you've used already.

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)