Skip to content

Instantly share code, notes, and snippets.

@euxoa
Created Apr 21, 2019
Embed
What would you like to do?
A simple Poisson state-space model with seasonality and overdispersion
data {
int<lower=1> n; // n = 2632
int y[n]; // number of shootings per day
}
transformed data {
int<lower=1,upper=365> yday[n];
for (i in 1:n) yday[i] = i % 365 + 1;
}
parameters {
vector[n] mu_innovations;
vector[365] seasonal_innovations; // yearly seasonality
vector[n] u_overdispersion;
real<lower=0> sigma_mu;
real<lower=0> sigma_yday;
real<lower=0> sigma_overdispersion;
real baseline;
}
transformed parameters {
// zero-mean seasonal term
vector[365] y_seasonal;
vector[n] mu;
{ vector[365] seasonal_with_trend;
real trend;
seasonal_with_trend = cumulative_sum(seasonal_innovations);
trend = seasonal_with_trend[365];
for (i in 1:365)
y_seasonal[i] = sigma_yday/100 * (seasonal_with_trend[i] - trend * i/365.0); }
mu = sigma_mu/100 * cumulative_sum(mu_innovations);
}
model {
seasonal_innovations ~ normal(0, 1);
mu_innovations ~ normal(0, 1);
u_overdispersion ~ normal(0, 1);
y ~ poisson_log(baseline + mu + y_seasonal[yday] + sigma_overdispersion * u_overdispersion);
sigma_mu ~ lognormal(-3.5 + log(100), 2);
sigma_yday ~ lognormal(-3.5 + log(100), 2);
sigma_overdispersion ~ lognormal(0, 2);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment