Hi all,
I’m using some code that was given by Michael Betancourt for the regularized horseshoe. The code runs fine, but I was wondering if there were any rules-of-thumb or general ways of thinking through how to choose the values in the \hl{transformed parameters} block.
Thanks,
David
modelString = "
data {
int <lower=1> n; // number of observations
int <lower=1> p; // number of predictors
real readscore[n]; // outcome
matrix[n,p] X; // inputs
}
transformed data {
real p0 = 5;
real slab_scale = 2.5;
real slab_scale2 = square(slab_scale);
real slab_df = 4;
real half_slab_df = 0.5 * slab_df;
}
parameters {
vector[p] beta_tilde;
vector<lower=0>[p] lambda;
real<lower=0> c2_tilde;
real<lower=0> tau_tilde;
real alpha;
real<lower=0> sigma;
}
transformed parameters {
vector [p] beta ; // regression coefficients
real tau0 = (p0 / (p - p0)) * (sigma / sqrt(1.0 * n));
real tau = tau0 * tau_tilde; // tau ~ cauchy(0, tau0)
// c2 ~ inv_gamma(half_slab_df, half_slab_df * slab_scale2)
// Implies that marginally beta ~ student_t(slab_df, 0, slab_scale)
real c2 = slab_scale2 * c2_tilde;
vector[p] lambda_tilde =
sqrt( c2 * square(lambda) ./ (c2 + square(tau) * square(lambda)) );
// beta ~ normal(0, tau * lambda_tilde)
beta = tau * lambda_tilde .* beta_tilde;
}
model {
beta_tilde ~ normal(0, 1);
lambda ~ cauchy(0, 1);
tau_tilde ~cauchy(0,1);
c2_tilde ~ inv_gamma(half_slab_df, half_slab_df);
alpha ~ normal(0, 2);
sigma ~ cauchy(0, 1);
readscore ~ normal(X * beta + alpha, sigma);
}
// For posterior predictive checking and loo cross-validation
generated quantities {
vector[n] readscore_rep;
vector[n] log_lik;
for (i in 1:n) {
readscore_rep[i] = normal_rng(alpha + X[i,:] * beta, sigma);
log_lik[i] = normal_lpdf(readscore[i] | alpha + X[i,:] * beta, sigma);
}
}
"