I am running this model and it seems that i always receive divergent transitions. Not just a few but seemingly all (about 2x my iterations). However, the parameter estimates are quite accurate and are pretty close to the true generating parameters. Is it safe to ignore this?
functions {
real log_Z_com_poisson_approx(real lambda, real nu) {
real log_mu = log(lambda^(1/nu));
real nu_mu = nu * exp(log_mu);
real nu2 = nu^2;
// first 4 terms of the residual series
real log_sum_resid = log1p(
nu_mu^(-1) * (nu2 - 1) / 24 +
nu_mu^(-2) * (nu2 - 1) / 1152 * (nu2 + 23) +
nu_mu^(-3) * (nu2 - 1) / 414720 * (5 * nu2^2 - 298 * nu2 + 11237)
);
return nu_mu + log_sum_resid -
((log(2 * pi()) + log_mu) * (nu - 1) / 2 + log(nu) / 2);
}
real compute_log_z(real lambda, real nu, real log_error) {
real z = negative_infinity();
real z_last = 0;
real j = 0;
while ((abs(z - z_last) > log_error) && j < 300) {
z_last = z;
z = log_sum_exp(z, j * log(lambda) - nu * lgamma(j+1));
j = j + 1;
if (j > 298) {
return log_Z_com_poisson_approx(lambda, nu);
}
}
//print(j);
return(z);
}
}
data {
int<lower=1> I; // number of observations
int<lower=1> J; // number of items
int<lower=1> N; // total number of counts
int<lower=0> Y[N]; // observed counts
int<lower=1, upper=I> person_id[N]; // person id for each observation
int<lower=1, upper=J> item_id[N]; // item id for each observation
}
parameters {
vector[I] ability;
vector[J] difficulty;
vector<lower=0>[J] nu;
}
transformed parameters {
real log_error = .001;
}
model {
// Priors
nu ~ lognormal(0, .25);
difficulty ~ normal(0,1);
ability ~ normal(0, 1);
// Likelihood
for (n in 1:N) {
real lambda = exp(nu[item_id[n]] * (ability[person_id[n]] + difficulty[item_id[n]]));
target += Y[n] * log(lambda) - nu[item_id[n]] * lgamma(Y[n]+1) - compute_log_z(lambda, nu[item_id[n]], log_error);
}
}