I am trying to fit multiple sets of parameters at once to a compartmental model. I was able to get a piece of code running that only fit one set of data to the model:
functions {
// This largely follows the deSolve package, but also includes the x_r and x_i variables.
// These variables are used in the background.
// params = [beta, sigma, gamma, delta]
// y = [S, E, I, R]
real[] SIER(real t,
real[] y,
real[] params,
real[] x_r,
int[] x_i) {
// Our compartments
real S = y[1];
real E = y[2];
real I = y[3];
real R = y[4];
// Our model parameters
real beta = params[1];
real sigma = params[2];
real gamma = params[3];
real delta = params[4];
real dydt[4];
dydt[1] = - beta * S * I;
dydt[2] = beta * S * I - sigma * E;
dydt[3] = sigma * E - (gamma + delta) * I;
dydt[4] = gamma * I;
return dydt;
}
}
data {
int<lower = 1> n_obs; // Number of days sampled
int<lower = 1> n_params; // Number of model parameters
int<lower = 1> n_difeq; // Number of differential equations in the system
int<lower = 100> n_pop; // This is to estimate the susceptible population
real<lower = 0> deaths[n_obs]; // The death data
real t0; // Initial time point (zero)
real ts[n_obs]; // Time index
}
transformed data {
real x_r[0];
int x_i[0];
}
parameters {
real<lower = 0, upper = 1> params[n_params]; // Model parameters
//int<lower = 0> S0; // Initial estimate of hosts susceptible
real<lower = 0> E0; // Initial estimate of infected hosts
real<lower = 0> I0; // Initial estimate of exposed hosts
real<lower = 0, upper = 1> prop_sus; // Estimated inverse of proportion of susceptible individuals
real<lower = 0> phi; // Variance of model
}
transformed parameters{
real y_hat[n_obs, n_difeq]; // Output from the ODE solver
real y0[n_difeq]; // Initial conditions for S, E, I, R
vector[n_obs] approx_deaths;
y0[1] = n_pop * prop_sus;
y0[2] = E0;
y0[3] = I0;
y0[4] = 0;
y_hat = integrate_ode_bdf(SIER, y0, t0, ts, params, x_r, x_i, 1e-9, 1e-5, 1e6);
approx_deaths = to_vector(y_hat[,3]) * params[4];
}
model {
params ~ beta(1, 20); //constrained to be 0-1
prop_sus ~ beta(10, 6); //constrained to be 0-1
E0 ~ normal(10, 2); //constrained to be positive.
I0 ~ normal(10, 2); //constrained to be positive.
phi ~ cauchy(15, 5); //constrained to be positive
deaths ~ normal(approx_deaths, phi);
}
generated quantities {
// Generate predicted data over the whole time series:
vector[n_obs] pred_deaths;
pred_deaths = to_vector(integrate_ode_bdf(SIER, y0, t0, ts, params, x_r, x_i, 1e-9, 1e-5, 1e6)[,3]) * params[4];
}
which works reasonably well. This is the code I am currently trying to run:
functions {
// This largely follows the deSolve package, but also includes the x_r and x_i variables.
// These variables are used in the background.
// params = [beta, sigma, gamma, delta]
// y = [S, E, I, R]
real[] SEIR(real t,
real[] y,
real[] params,
real[] x_r,
int[] x_i) {
// Our compartments
real S = y[1];
real E = y[2];
real I = y[3];
real R = y[4];
// Our model parameters
real beta = params[1];
real sigma = params[2];
real gamma = params[3];
real delta = params[4];
real dydt[4];
dydt[1] = - beta * S * I;
dydt[2] = beta * S * I - sigma * E;
dydt[3] = sigma * E - (gamma + delta) * I;
dydt[4] = gamma * I;
return dydt;
}
}
data {
int<lower = 1> n_obs; // Number of days sampled
int<lower = 1> n_states; // Number of state data provided
int<lower = 1> n_params; // Number of model parameters
int<lower = 1> n_difeq; // Number of differential equations in the system
int<lower = 100> n_pop[n_states]; // This is to estimate the susceptible population
real<lower = 0> deaths[n_states, n_obs]; // The death data
real t0; // Initial time point (zero)
real ts[n_obs]; // Time index
}
transformed data {
real x_r[0];
int x_i[0];
}
parameters {
real<lower = 0, upper = 1> params[n_states, n_params]; // Model parameters
//int<lower = 0> S0; // Initial estimate of hosts susceptible
real<lower = 0> I0[n_states]; // Initial estimate of infected hosts
real<lower = 0> E0[n_states]; // Initial estimate of exposed hosts
real<lower = 0, upper = 1> prop_sus[n_states]; // Estimated inverse of proportion of susceptible individuals
real<lower = 0> phi[n_states];
}
transformed parameters{
real y_hat[n_states, n_obs, n_difeq]; // Output from the ODE solver
real y0[n_states, n_difeq]; // Initial conditions for S, E, I, R
vector[n_obs] approx_deaths[n_states];
for (i in 1:n_states) {
y0[i][1] = n_pop[i] * prop_sus[i];
y0[i][2] = E0[i];
y0[i][3] = I0[i];
y0[i][4] = 0;
y_hat[i] = integrate_ode_bdf(SEIR, y0[i], t0, ts, params[i], x_r, x_i, 1e-9, 1e-5, 1e6);
approx_deaths[i] = to_vector(y_hat[i][,3]) * params[i][4];
}
}
model {
prop_sus ~ beta(10, 6); //constrained to be 0-1
E0 ~ normal(10, 2); //constrained to be positive.
I0 ~ normal(10, 2); //constrained to be positive.
phi ~ cauchy(15, 5); //constrained to be positive
for (i in 1:n_states) {
params[i] ~ beta(1, 20); //constrained to be 0-1
deaths[i] ~ normal(approx_deaths[i], phi[i]);
}
}
generated quantities {
// Generate predicted data over the whole time series:
vector[n_obs] pred_deaths[n_states];
for (i in 1:n_states) {
pred_deaths[i] = to_vector(integrate_ode_bdf(SEIR, y0[i], t0, ts, params[i], x_r, x_i, 1e-9, 1e-5, 1e6)[,3]) * params[i][4];
}
}
Thanks!