I just upgraded to CmdStan 2.3 and when I try to sample the prior model, I get lots of NA values when I inspect “pred” (generated quantities). This did not happen when I run the same model with 2.29.
I’m not sure if this is a bug or if the new version is just finding a problem with the model that 2.29 was unable to detect.
I’m attaching the model and the R workspace with the data
Model:
data {
//Dimensions
int<lower=1> N;
int<lower=1> N_site;
int<lower=1> N_plot;
int<lower=1> N_cover;
//Response variable
array[N] int<lower=0> germ;
//Explanatory variables
vector[N] SM;
vector[N] ST;
array[N] int<lower=1,upper=N_cover> cover_id;
//group variables
array[N] int<lower=1> site_id;
array[N] int<lower=1> plot_id;
}
parameters {
//intercept
real a;
//varying intercepts
vector[N_site] a_site_nc;
real<lower=0> sigma_site;
vector[N_plot] a_plot;
vector[N_cover] b_sm_nc;
real mu_sm;
real<lower=0>sigma_sm;
vector[N_cover] b_st_nc;
real mu_st;
real<lower=0>sigma_st;
//BB
real<lower=0> phi;
}
transformed parameters{
real theta;
vector[N_site] a_site = 0 + sigma_site*a_site_nc;
vector[N_cover] b_st = mu_st+sigma_st*b_st_nc;
vector[N_cover] b_sm = mu_sm+sigma_sm*b_sm_nc;
vector[N] p = inv_logit(a + a_site[site_id]+a_plot[plot_id]+
b_sm[cover_id].*SM+b_st[cover_id].*ST);
theta = phi + 2;
}
model {
//Likelihood
//germ~beta_binomial(6, p * theta, (1-p) * theta);
b_sm_nc~normal(0,1);
mu_sm~normal(0,0.1);
sigma_sm~exponential(3);
b_st_nc~normal(0,1);
mu_st~normal(0,0.1);
sigma_st~exponential(3);
//priors for intercepts
a~normal(0,0.1);
a_site_nc~normal(0,1);
sigma_site~exponential(3);
a_plot~normal(0,0.1);
//bb
phi ~ exponential(1);
}
generated quantities {
vector[N] pred;
for (i in 1:N) {
pred[i] = beta_binomial_rng(6, p[i]*theta, (1-p[i])*theta);
}
}
Runn the code
library(cmdstanr)
mod<-cmdstan_model("prior_model.stan")
fit_prior<-mod$sample(data=data_list,
chains=4,
parallel_chains = 4)
prior_samples<-fit_prior$draws(c("pred"),format="matrix")
workspace.RData (3.6 KB)
prior_model.stan (1.5 KB)