Hello!
I have a hierarchical model taking approximately 3 minutes to fit.
//
// This Stan program defines a hierarchical model, with a
// vector of values isotopic 'y' modeled as normally distributed
// with hierarchical slopes between C/N and y
// being dependant upon soil C content
//
data {
// Integers
int<lower=0,upper=1> prior_only;
int<lower=0> N; //nb of data points
int<lower=0> U; //nb of units
// Response
vector[N] y;
// Covariates
vector[U] SOC;
vector[N] CN;
array[N] int unit;
// Priors
real lam_sigma_beta_1; real lam_sigma_beta_2;
real lam_sigma;
real sd_gamma_0;
real sd_gamma;
real mean_beta_0; real sd_beta_0;
}
parameters {
// Hyperparameters
cholesky_factor_corr[2] rho; //correlation among hierarchical effects
vector<lower=0>[2] sigma_beta; //sd of hierarchical effects
real<multiplier=sd_gamma_0> gamma_0; //int of the slope-SOC relationship
real<multiplier=sd_gamma> gamma; //slope of the slope-SOC relationship
// Parameters
real<offset=mean_beta_0,multiplier=sd_beta_0> beta_0; //int of the isotope-CN relationship
matrix[U,2] beta_raw; //slope of the isotope-CN relationship
real<lower=0> sigma; //residual sd
}
transformed parameters{
// Hyperparameters
matrix[U,2] mu_beta = append_col(rep_vector(0,U), gamma_0 + SOC * gamma); //linear predictor of slopes
matrix[2,2] Sigma_beta = diag_post_multiply(rho, sigma_beta); //cov matrix of hierarchical effects
// Parameters
matrix[U,2] beta = mu_beta + beta_raw * Sigma_beta; //final hierarchical parameters
// Linear predictor
vector[N] mu = beta_0 + beta[unit,1] + CN .* beta[unit,2];
}
model {
// Hyperpriors
rho ~ lkj_corr_cholesky(2);
sigma_beta[1] ~ exponential(lam_sigma_beta_1);
sigma_beta[2] ~ exponential(lam_sigma_beta_2);
gamma_0 ~ normal(0, sd_gamma_0);
gamma ~ normal(0, sd_gamma);
// Priors
beta_0 ~ normal(mean_beta_0, sd_beta_0);
to_vector(beta_raw) ~ std_normal();
sigma ~ exponential(lam_sigma);
// Likelihood
if(prior_only == 0){
y ~ normal(mu, sigma);
}
}
generated quantities{
real y_rep[N] = normal_rng(mu, sigma);
}
When I had ubuntu OS, this very simple standalone generated quantities took some seconds to compute.
//
// This Stan program defines the generated quantities for a hierarchical model, with a
// vector of values isotopic 'y' modeled as normally distributed
// with hierarchical slopes between C/N and y
// being dependant upon soil C content
//
data {
// Integers
int<lower=0,upper=1> prior_only;
int<lower=0> N; //nb of data points
int<lower=0> N_pred; //nb of C/N values for predictions
int<lower=0> U; //nb of units
int<lower=0> U_pred; //nb of SOC values for predictions
// Response
vector[N] y;
// Covariates
vector[U] SOC;
vector[U_pred] SOC_pred;
vector[N] CN;
vector[N_pred] CN_pred;
int unit[N];
int unit_pred[N_pred];
// Priors
real lam_sigma_beta_1; real lam_sigma_beta_2;
real lam_sigma;
real sd_gamma_0;
real sd_gamma;
real mean_beta_0; real sd_beta_0;
}
parameters {
// Hyperparameters
cholesky_factor_corr[2] rho; //correlation among hierarchical effects
vector<lower=0>[2] sigma_beta; //sd of hierarchical effects
real<multiplier=sd_gamma_0> gamma_0; //int of the slope-SOC relationship
real<multiplier=sd_gamma> gamma; //slope of the slope-SOC relationship
// Parameters
real<offset=mean_beta_0,multiplier=sd_beta_0> beta_0; //int of the isotope-CN relationship
matrix[U,2] beta_raw; //slope of the isotope-CN relationship
real<lower=0> sigma; //residual sd
}
transformed parameters{
// Hyperparameters
matrix[U,2] mu_beta = append_col(rep_vector(0,U), gamma_0 + SOC * gamma); //linear predictor of slopes
matrix[2,2] Sigma_beta = diag_post_multiply(rho, sigma_beta); //cov matrix of hierarchical effects
// Parameters
matrix[U,2] beta = mu_beta + beta_raw * Sigma_beta; //final hierarchical parameters
// Linear predictor
vector[N] mu = beta_0 + beta[unit,1] + CN .* beta[unit,2];
}
generated quantities{
vector[U_pred] mu_beta_pred = gamma_0 + SOC_pred * gamma;
vector[N_pred] mu_pred = beta_0 + beta[unit_pred,1] + CN_pred .* beta[unit_pred,2];
}
Now that I switched to windows, it takes almost the same time as model fitting! I use as many cores in the GQ calls as in the fitting call.
Does anyone already noticed that? Is there a way to make GQ faster?
I am using cmstanR 0.5.2 with the last cmdstan.
Thank you very much!
All the best,
Lucas