Hello,
I’m very new to Stan, but am excited to get started. I am trying to convert a random effects model from pymc3 into Stan.
I have been able to use pystan to compile a random effects model and draw posterior samples from the fitted distribution. See below for the code:
// hierarchical model with non-centered parameterization of mu
data {
int<lower=1> N; //Number of observations
int<lower=0> J; //Number of groups
int<lower=0> K; //Number of group predictors
int<lower=0> O; //Number of observation predictors
int<lower=0, upper=J> group_idx[N]; // Maps each observation to each group
real response[N]; // The outcome of the data generation process
matrix[N,O] observation_data; //covariate data of the observations
matrix[J,K] group_data; //covariate data of the groups
}
parameters {
real<lower=2> nu;
real<lower=0> inv_phi;
vector[O] beta;
vector[K] gamma;
vector[J] group_mu_raw;
real<lower=0> sigma_group_mu;
}
transformed parameters {
vector[J] group_mu = group_data * gamma + sigma_group_mu * group_mu_raw;
vector[N] mu =
group_mu[group_idx] +
observation_data * beta;
real phi = 1 / inv_phi;
}
model {
response ~ student_t(nu, mu, phi);
nu ~ gamma(2, 0.1);
inv_phi ~ student_t(4, 0, 5);
beta ~ student_t(4, 0, 5);
gamma ~ student_t(4, 0, 5);
group_mu_raw ~ student_t(4, 0, 5); // std_normal()
sigma_group_mu ~ student_t(4, 0, 5);
}
generated quantities {
real y_rep[N];
for (n in 1:N) {
y_rep[n] = student_t_rng(nu, mu[n], phi);
}
}
However, I now want to generate predictions from out-of-sample data. To do this, I have created a second Stan model that only contains the data and generated quantities block as recommended elsewhere. See below for code:
// hierarchical model with non-centered parameterization of mu
data {
int<lower=1> N; //Number of observations
int<lower=0> J; //Number of groups
int<lower=0> K; //Number of group predictors
int<lower=0> O; //Number of observation predictor
int<lower=0, upper=J> group_idx[N]; // Maps each observation to each group
real response[N]; // The outcome of the data generation process
matrix[N,O] observation_data; //covariate data of the observations
matrix[J,K] group_data; //covariate data of the groups
int<lower=0> max_group_idx; //Highest
int N_samples; // number of samples drawn from posterior distributions
real nu[N_samples]; // nu parameter from student t likelihood
real phi[N_samples]; // standard deviation parameter from student t likelihood
matrix[N_samples, O] beta; // covariate matrix draws for observations
matrix[N_samples, K] gamma; // covariate matrix draws for groups
matrix[N_samples, max_group_idx] group_mu; //posterior draws of the group mu parameter
real sigma_group_mu[N_samples]; // posterior draws of the hierarchal variance parameter.
}
parameters {
}
model {
}
generated quantities {
real y_pred[N, N_samples];
row_vector[N_samples] group_mu_rep;
matrix[N, N_samples] mu;
vector[O] observation_data_vector;
for(n in 1:N) {
if (group_idx[n] > max_group_idx) {
for (s in 1:N_samples) {
group_mu_rep[s] = gamma[s] * to_vector(group_data[group_idx[n]]) + sigma_group_mu[s] * student_t_rng(4, 0, 1);
}
}
else {
group_mu_rep = to_row_vector(group_mu[,group_idx[n]]);
}
mu[n] = group_mu_rep + to_row_vector(beta * to_vector(observation_data[n]));
y_pred[n] = student_t_rng(nu, mu[n], phi);
}
}
The problem now is that it takes an extremely long time to compile the model and draw samples from the predicted distribution.
It also feels kind of a hack to specify iter=1 when sampling.
pred_model.sampling(data=model_data, chains=1, iter=1, algorithm="Fixed_param")
Are there any recommendations for specifying the generated_quantities block in the prediction or how to speed this operation up?