Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
data {
int<lower=1> N_ages;
vector[N_ages] ages;
int K[N_ages];
int N[N_ages];
}
transformed data {
vector[N_ages] mu;
mu <- rep_vector(0, N_ages);
}
parameters {
real<lower=0> theta_1;
real<lower=0> theta_2;
real<lower=0> theta_3;
real<lower=0> theta_4;
vector[N_ages] y;
}
model {
matrix[N_ages, N_ages] Sigma;
for (i in 1:N_ages)
for (j in i:N_ages) {
Sigma[i, j] <- theta_1 * exp(- theta_2 * square(ages[i] - ages[j])) + theta_3 + theta_4 * ages[i] * ages[j];
}
for (i in 1:N_ages)
for (j in (i+1):N_ages)
Sigma[j, i] <- Sigma[i, j];
// Priors
theta_1 ~ cauchy(0, 5);
theta_2 ~ cauchy(0, 5);
theta_3 ~ cauchy(0, 5);
theta_4 ~ cauchy(0, 5);
y ~ multi_normal(mu, Sigma);
// Likelihood
K ~ binomial_logit(N, y);
}
generated quantities {
vector[N_ages] p_post;
vector[N_ages] kdn_post;
for (i in 1:N_ages) {
p_post[i] <- inv_logit(y[i]);
kdn_post[i] <- 1.0 * binomial_rng(N[i], p_post[i]) / N[i];
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment