Howdy! I’ve been trying to code up a basic structural time series (state-space) model following the code in this tutorial but modifying it for a negative binomial likelihood.
Unfortunately, when I run the model it always seems to have less than 4 chains run - sometimes 1, sometimes 3, whatever. I have tried to specify chains=1 to debug, but every time I have tried that, it runs, so I can’t find the error.
I have included a reproducible example below. (Note that the neg_binomial_2_log_safe_rng is from Numerical stability of GPs with negative binomial likelihood - #5 by jonah because I was having the same errors as mentioned in that thread.)
I’m hoping someone can point out what is wrong with my model. I’ve probably missed something obvious in my adaption from the Gaussian model to negative binomial. Thanks!
library(rstan)
library(bayesplot)
set.seed(1435)
T <- 100
y <- rep(NA, T)
shape <- 20
mu_err <- rnorm(T, 0, 0.035)
delta_err <- rnorm(T, 0, 0.025)
mu <- rep(NA, T)
delta <- rep(NA, T)
mu[1] <- mu_err[1]
delta[1] <- delta_err[1]
for (t in 2:T) {
mu[t] <- mu[t-1] + delta[t-1] + mu_err[t]
delta[t] <- delta[t-1] + delta_err[t]
}
y <- rnbinom(n=T, mu=exp(mu), size=shape)
time <- 1:T
plot(y ~ time)
#stan data
d1 <- list(T = T, #total number of observations, this is also the total time in the time series
y = y #response variable
)
#stan code
stan_code1 <- "
functions {
int neg_binomial_2_log_safe_rng(real eta, real phi) {
real gamma_rate = gamma_rng(phi, phi / exp(eta));
if (gamma_rate > exp(20.7)) gamma_rate = exp(20.7); // Jonah Gabry thinks this is the max value before overflow but hasn't double checked
return poisson_rng(gamma_rate);
}
}
data {
int<lower=0> T;
int y[T];
}
parameters {
vector[T] mu_err;
vector[T] delta_err;
real<lower=0> s_slope;
real<lower=0> s_level;
real<lower=0> shape;
}
transformed parameters {
vector[T] mu;
vector[T] delta;
mu[1] = mu_err[1];
delta[1] = delta_err[1];
for (t in 2:T) {
mu[t] = mu[t-1] + delta[t-1] + s_level * mu_err[t];
delta[t] = delta[t-1] + s_slope * delta_err[t];
}
}
model {
for (t in 1:T) {
target += neg_binomial_2_log_lpmf(y[t] | mu[t], shape);
}
target += gamma_lpdf(shape | 0.01, 0.01);
target += normal_lpdf(s_slope | 0, 0.5);
target += normal_lpdf(s_level | 0, 0.5);
target += std_normal_lpdf(mu_err);
target += std_normal_lpdf(delta_err);
}
generated quantities {
vector[T] preds;
for (t in 1:T) {
preds[t] = neg_binomial_2_log_safe_rng(mu[t], shape);
}
}
"
#fit model
fit1 <- stan(model_code = stan_code1, data = d1,
chains = 4, iter = 2000, warmup = 1000,
thin = 1, cores = 4)
print(fit1)
ypreds <- extract(fit1)$preds
ppc_bars(y, yrep = ypreds[1:30, ])
#fit same model with one chain (and less iter for speed sake)
fit2 <- stan(model_code = stan_code1, data = d1,
chains = 1, iter = 1000, warmup = 500,
thin = 1, cores = 1)
