I share my RStan code. I want to implement a 10 age-structured SIR model.
I don’t know anything about the priors.
However just to simulate which trajectories I can obtain (avoiding the likelihood) the gradient evaluation takes about 0.4 seconds.
Do you have any suggestions?
functions{
real[] sir(real t, // the is a real variable
real[] y, // the vector field is made with real funtions
real[] theta, // the parameters vector
real[] x_r, // these are vabiables that depend from the parameters
int[] x_i
){ // integer variables
real S_0_9 = y[1];
real I_0_9 = y[2];
real R_0_9 = y[3];
real N_0_9 = x_i[1];
real S_10_19 = y[4];
real I_10_19 = y[5];
real R_10_19 = y[6];
real N_10_19 = x_i[2];
real S_20_29 = y[7];
real I_20_29 = y[8];
real R_20_29 = y[9];
real N_20_29 = x_i[3];
real S_30_39 = y[10];
real I_30_39 = y[11];
real R_30_39 = y[12];
real N_30_39 = x_i[4];
real S_40_49 = y[13];
real I_40_49 = y[14];
real R_40_49 = y[15];
real N_40_49 = x_i[5];
real S_50_59 = y[16];
real I_50_59 = y[17];
real R_50_59 = y[18];
real N_50_59 = x_i[6];
real S_60_69 = y[19];
real I_60_69 = y[20];
real R_60_69 = y[21];
real N_60_69 = x_i[7];
real S_70_79 = y[22];
real I_70_79 = y[23];
real R_70_79 = y[24];
real N_70_79 = x_i[8];
real S_80_89 = y[25];
real I_80_89 = y[26];
real R_80_89 = y[27];
real N_80_89 = x_i[9];
real S_over90 = y[28];
real I_over90 = y[29];
real R_over90 = y[30];
real N_over90 = x_i[10];
real beta = theta[1];
real gamma_0_9 = theta[2];
real gamma_10_19 = theta[3];
real gamma_20_29 = theta[4];
real gamma_30_39 = theta[5];
real gamma_40_49 = theta[6];
real gamma_50_59 = theta[7];
real gamma_60_69 = theta[8];
real gamma_70_79 = theta[9];
real gamma_80_89 = theta[10];
real gamma_over90 = theta[11];
real c_0_9_0_9 = x_r[1];
real c_0_9_10_19 = x_r[2];
real c_0_9_20_29 = x_r[3];
real c_0_9_30_39 = x_r[4];
real c_0_9_40_49 = x_r[5];
real c_0_9_50_59 = x_r[6];
real c_0_9_60_69 = x_r[7];
real c_0_9_70_79 = x_r[8];
real c_0_9_80_89 = x_r[9];
real c_0_9_over90 = x_r[10];
real c_10_19_0_9 = x_r[11];
real c_10_19_10_19 = x_r[12];
real c_10_19_20_29 = x_r[13];
real c_10_19_30_39 = x_r[14];
real c_10_19_40_49 = x_r[15];
real c_10_19_50_59 = x_r[16];
real c_10_19_60_69 = x_r[17];
real c_10_19_70_79 = x_r[18];
real c_10_19_80_89 = x_r[19];
real c_10_19_over90 = x_r[20];
real c_20_29_0_9 = x_r[21];
real c_20_29_10_19 = x_r[22];
real c_20_29_20_29 = x_r[23];
real c_20_29_30_39 = x_r[24];
real c_20_29_40_49 = x_r[25];
real c_20_29_50_59 = x_r[26];
real c_20_29_60_69 = x_r[27];
real c_20_29_70_79 = x_r[28];
real c_20_29_80_89 = x_r[20];
real c_20_29_over90 = x_r[30];
real c_30_39_0_9 = x_r[31];
real c_30_39_10_19 = x_r[32];
real c_30_39_20_29 = x_r[33];
real c_30_39_30_39 = x_r[34];
real c_30_39_40_49 = x_r[35];
real c_30_39_50_59 = x_r[36];
real c_30_39_60_69 = x_r[37];
real c_30_39_70_79 = x_r[38];
real c_30_39_80_89 = x_r[39];
real c_30_39_over90 = x_r[40];
real c_40_49_0_9 = x_r[41];
real c_40_49_10_19 = x_r[42];
real c_40_49_20_29 = x_r[43];
real c_40_49_30_39 = x_r[44];
real c_40_49_40_49 = x_r[45];
real c_40_49_50_59 = x_r[46];
real c_40_49_60_69 = x_r[47];
real c_40_49_70_79 = x_r[48];
real c_40_49_80_89 = x_r[49];
real c_40_49_over90 = x_r[50];
real c_50_59_0_9 = x_r[51];
real c_50_59_10_19 = x_r[52];
real c_50_59_20_29 = x_r[53];
real c_50_59_30_39 = x_r[54];
real c_50_59_40_49 = x_r[55];
real c_50_59_50_59 = x_r[56];
real c_50_59_60_69 = x_r[57];
real c_50_59_70_79 = x_r[58];
real c_50_59_80_89 = x_r[59];
real c_50_59_over90 = x_r[60];
real c_60_69_0_9 = x_r[61];
real c_60_69_10_19 = x_r[62];
real c_60_69_20_29 = x_r[63];
real c_60_69_30_39 = x_r[64];
real c_60_69_40_49 = x_r[65];
real c_60_69_50_59 = x_r[66];
real c_60_69_60_69 = x_r[67];
real c_60_69_70_79 = x_r[68];
real c_60_69_80_89 = x_r[69];
real c_60_69_over90 = x_r[70];
real c_70_79_0_9 = x_r[71];
real c_70_79_10_19 = x_r[72];
real c_70_79_20_29 = x_r[73];
real c_70_79_30_39 = x_r[74];
real c_70_79_40_49 = x_r[75];
real c_70_79_50_59 = x_r[76];
real c_70_79_60_69 = x_r[77];
real c_70_79_70_79 = x_r[78];
real c_70_79_80_89 = x_r[79];
real c_70_79_over90 = x_r[80];
real c_80_89_0_9 = x_r[81];
real c_80_89_10_19 = x_r[82];
real c_80_89_20_29 = x_r[83];
real c_80_89_30_39 = x_r[84];
real c_80_89_40_49 = x_r[85];
real c_80_89_50_59 = x_r[86];
real c_80_89_60_69 = x_r[87];
real c_80_89_70_79 = x_r[88];
real c_80_89_80_89 = x_r[89];
real c_80_89_over90 = x_r[90];
real c_over90_0_9 = x_r[91];
real c_over90_10_19 = x_r[92];
real c_over90_20_29 = x_r[93];
real c_over90_30_39 = x_r[94];
real c_over90_40_49 = x_r[95];
real c_over90_50_59 = x_r[96];
real c_over90_60_69 = x_r[97];
real c_over90_70_79 = x_r[98];
real c_over90_80_89 = x_r[99];
real c_over90_over90 = x_r[100];
real lambda_0_9 = beta * (c_0_9_0_9 * I_0_9/N_0_9 + c_0_9_10_19 * I_10_19/N_10_19 + c_0_9_20_29 * I_20_29/N_20_29
+ c_0_9_30_39 * I_30_39/N_30_39 + c_0_9_40_49 * I_40_49/N_40_49 + c_0_9_50_59 * I_50_59/N_50_59
+ c_0_9_60_69 * I_60_69/N_60_69 + c_0_9_70_79 * I_70_79/N_70_79 + c_0_9_80_89 * I_80_89/N_80_89
+ c_0_9_over90 * I_over90/N_over90) ;
real lambda_10_19 = beta * (c_10_19_0_9 * I_0_9/N_0_9 + c_10_19_10_19 * I_10_19/N_10_19 + c_10_19_20_29 * I_20_29/N_20_29
+ c_10_19_30_39 * I_30_39/N_30_39 + c_10_19_40_49 * I_40_49/N_40_49 + c_10_19_50_59 * I_50_59/N_50_59
+ c_10_19_60_69 * I_60_69/N_60_69 + c_10_19_70_79 * I_70_79/N_70_79 + c_10_19_80_89 * I_80_89/N_80_89
+ c_10_19_over90 * I_over90/N_over90) ;
real lambda_20_29 = beta * (c_20_29_0_9 * I_0_9/N_0_9 + c_20_29_10_19 * I_10_19/N_10_19 + c_20_29_20_29 * I_20_29/N_20_29
+ c_20_29_30_39 * I_30_39/N_30_39 + c_20_29_40_49 * I_40_49/N_40_49 + c_20_29_50_59 * I_50_59/N_50_59
+ c_20_29_60_69 * I_60_69/N_60_69 + c_20_29_70_79 * I_70_79/N_70_79 + c_20_29_80_89 * I_80_89/N_80_89
+ c_20_29_over90 * I_over90/N_over90) ;
real lambda_30_39 = beta * (c_30_39_0_9 * I_0_9/N_0_9 + c_30_39_10_19 * I_10_19/N_10_19 + c_30_39_20_29 * I_20_29/N_20_29
+ c_30_39_30_39 * I_30_39/N_30_39 + c_30_39_40_49 * I_40_49/N_40_49 + c_30_39_50_59 * I_50_59/N_50_59
+ c_30_39_60_69 * I_60_69/N_60_69 + c_30_39_70_79 * I_70_79/N_70_79 + c_30_39_80_89 * I_80_89/N_80_89
+ c_30_39_over90 * I_over90/N_over90) ;
real lambda_40_49 = beta * (c_40_49_0_9 * I_0_9/N_0_9 + c_40_49_10_19 * I_10_19/N_10_19 + c_40_49_20_29 * I_20_29/N_20_29
+ c_40_49_30_39 * I_30_39/N_30_39 + c_40_49_40_49 * I_40_49/N_40_49 + c_40_49_50_59 * I_50_59/N_50_59
+ c_40_49_60_69 * I_60_69/N_60_69 + c_40_49_70_79 * I_70_79/N_70_79 + c_40_49_80_89 * I_80_89/N_80_89
+ c_40_49_over90 * I_over90/N_over90) ;
real lambda_50_59 = beta * (c_50_59_0_9 * I_0_9/N_0_9 + c_50_59_10_19 * I_10_19/N_10_19 + c_50_59_20_29 * I_20_29/N_20_29
+ c_50_59_30_39 * I_30_39/N_30_39 + c_50_59_40_49 * I_40_49/N_40_49 + c_50_59_50_59 * I_50_59/N_50_59
+ c_50_59_60_69 * I_60_69/N_60_69 + c_50_59_70_79 * I_70_79/N_70_79 + c_50_59_80_89 * I_80_89/N_80_89
+ c_50_59_over90 * I_over90/N_over90) ;
real lambda_60_69 = beta * (c_60_69_0_9 * I_0_9/N_0_9 + c_60_69_10_19 * I_10_19/N_10_19 + c_60_69_20_29 * I_20_29/N_20_29
+ c_60_69_30_39 * I_30_39/N_30_39 + c_60_69_40_49 * I_40_49/N_40_49 + c_60_69_50_59 * I_50_59/N_50_59
+ c_60_69_60_69 * I_60_69/N_60_69 + c_60_69_70_79 * I_70_79/N_70_79 + c_60_69_80_89 * I_80_89/N_80_89
+ c_60_69_over90 * I_over90/N_over90) ;
real lambda_70_79 = beta * (c_70_79_0_9 * I_0_9/N_0_9 + c_70_79_10_19 * I_10_19/N_10_19 + c_70_79_20_29 * I_20_29/N_20_29
+ c_70_79_30_39 * I_30_39/N_30_39 + c_70_79_40_49 * I_40_49/N_40_49 + c_70_79_50_59 * I_50_59/N_50_59
+ c_70_79_60_69 * I_60_69/N_60_69 + c_70_79_70_79 * I_70_79/N_70_79 + c_70_79_80_89 * I_80_89/N_80_89
+ c_70_79_over90 * I_over90/N_over90) ;
real lambda_80_89 = beta * (c_80_89_0_9 * I_0_9/N_0_9 + c_80_89_10_19 * I_10_19/N_10_19 + c_80_89_20_29 * I_20_29/N_20_29
+ c_80_89_30_39 * I_30_39/N_30_39 + c_80_89_40_49 * I_40_49/N_40_49 + c_80_89_50_59 * I_50_59/N_50_59
+ c_80_89_60_69 * I_60_69/N_60_69 + c_80_89_70_79 * I_70_79/N_70_79 + c_80_89_80_89 * I_80_89/N_80_89
+ c_80_89_over90 * I_over90/N_over90) ;
real lambda_over90 = beta * (c_over90_0_9 * I_0_9/N_0_9 + c_over90_10_19 * I_10_19/N_10_19 + c_over90_20_29 * I_20_29/N_20_29
+ c_over90_30_39 * I_30_39/N_30_39 + c_over90_40_49 * I_40_49/N_40_49 + c_over90_50_59 * I_50_59/N_50_59
+ c_over90_60_69 * I_60_69/N_60_69 + c_over90_70_79 * I_70_79/N_70_79 + c_over90_80_89 * I_80_89/N_80_89
+ c_over90_over90 * I_over90/N_over90) ;
real dS_0_9_dt = - lambda_0_9 * S_0_9 ;
real dI_0_9_dt = + lambda_0_9 * S_0_9 - gamma_0_9 * I_0_9;
real dR_0_9_dt = + gamma_0_9 * I_0_9;
real dS_10_19_dt = - lambda_10_19 * S_10_19 ;
real dI_10_19_dt = + lambda_10_19 * S_10_19 - gamma_10_19 * I_10_19;
real dR_10_19_dt = + gamma_10_19 * I_10_19;
real dS_20_29_dt = - lambda_20_29 * S_20_29 ;
real dI_20_29_dt = + lambda_20_29 * S_20_29 - gamma_20_29 * I_20_29;
real dR_20_29_dt = + gamma_20_29 * I_20_29;
real dS_30_39_dt = - lambda_30_39 * S_30_39 ;
real dI_30_39_dt = + lambda_30_39 * S_30_39 - gamma_30_39 * I_30_39;
real dR_30_39_dt = + gamma_30_39 * I_30_39;
real dS_40_49_dt = - lambda_40_49 * S_40_49 ;
real dI_40_49_dt = + lambda_40_49 * S_40_49 - gamma_40_49 * I_40_49;
real dR_40_49_dt = + gamma_40_49 * I_40_49;
real dS_50_59_dt = - lambda_50_59 * S_50_59 ;
real dI_50_59_dt = + lambda_50_59 * S_50_59 - gamma_50_59 * I_50_59;
real dR_50_59_dt = + gamma_50_59 * I_50_59;
real dS_60_69_dt = - lambda_60_69 * S_60_69 ;
real dI_60_69_dt = + lambda_60_69 * S_60_69 - gamma_60_69 * I_60_69;
real dR_60_69_dt = + gamma_60_69 * I_60_69;
real dS_70_79_dt = - lambda_70_79 * S_70_79 ;
real dI_70_79_dt = + lambda_70_79 * S_70_79 - gamma_70_79 * I_70_79;
real dR_70_79_dt = + gamma_70_79 * I_70_79;
real dS_80_89_dt = - lambda_80_89 * S_80_89 ;
real dI_80_89_dt = + lambda_80_89 * S_80_89 - gamma_80_89 * I_80_89;
real dR_80_89_dt = + gamma_80_89 * I_80_89;
real dS_over90_dt = - lambda_over90 * S_over90 ;
real dI_over90_dt = + lambda_over90 * S_over90 - gamma_over90 * I_over90;
real dR_over90_dt = + gamma_over90 * I_over90;
return{dS_0_9_dt, dI_0_9_dt, dR_0_9_dt, dS_10_19_dt, dI_10_19_dt, dR_10_19_dt,
dS_20_29_dt, dI_20_29_dt, dR_20_29_dt,dS_30_39_dt, dI_30_39_dt, dR_30_39_dt,
dS_40_49_dt, dI_40_49_dt, dR_40_49_dt,dS_50_59_dt, dI_50_59_dt, dR_50_59_dt,
dS_60_69_dt, dI_60_69_dt, dR_60_69_dt,dS_70_79_dt, dI_70_79_dt, dR_70_79_dt,
dS_80_89_dt, dI_80_89_dt, dR_80_89_dt,dS_over90_dt, dI_over90_dt, dR_over90_dt};
}
}
// -----------------------------------------------------------------------------------------
data { int<lower = 1> n_days;
int compute_likelihood;
real y0[30];
real initial_time;
real ts[n_days];
int N_data[10];
int I_0_9_data[n_days];
int I_10_19_data[n_days];
int I_20_29_data[n_days];
int I_30_39_data[n_days];
int I_40_49_data[n_days];
int I_50_59_data[n_days];
int I_60_69_data[n_days];
int I_70_79_data[n_days];
int I_80_89_data[n_days];
int I_over90_data[n_days];
real contact_matrix[10,10];
}
// ---------------------------------------------------------------------------------
transformed data {
real x_r[100] = {contact_matrix[1,1],contact_matrix[1,2],contact_matrix[1,3],contact_matrix[1,4],contact_matrix[1,5],
contact_matrix[1,6],contact_matrix[1,7],contact_matrix[1,8],contact_matrix[1,9],contact_matrix[1,10],
contact_matrix[2,1],contact_matrix[2,2],contact_matrix[2,3],contact_matrix[2,4],contact_matrix[2,5],
contact_matrix[2,6],contact_matrix[2,7],contact_matrix[2,8],contact_matrix[2,9],contact_matrix[2,10],
contact_matrix[3,1],contact_matrix[3,2],contact_matrix[3,3],contact_matrix[3,4],contact_matrix[3,5],
contact_matrix[3,6],contact_matrix[3,7],contact_matrix[3,8],contact_matrix[3,9],contact_matrix[3,10],
contact_matrix[4,1],contact_matrix[4,2],contact_matrix[4,3],contact_matrix[4,4],contact_matrix[4,5],
contact_matrix[4,6],contact_matrix[4,7],contact_matrix[4,8],contact_matrix[4,9],contact_matrix[4,10],
contact_matrix[5,1],contact_matrix[5,2],contact_matrix[5,3],contact_matrix[5,4],contact_matrix[5,5],
contact_matrix[5,6],contact_matrix[5,7],contact_matrix[5,8],contact_matrix[5,9],contact_matrix[5,10],
contact_matrix[6,1],contact_matrix[6,2],contact_matrix[6,3],contact_matrix[6,4],contact_matrix[6,5],
contact_matrix[6,6],contact_matrix[6,7],contact_matrix[6,8],contact_matrix[6,9],contact_matrix[6,10],
contact_matrix[7,1],contact_matrix[7,2],contact_matrix[7,3],contact_matrix[7,4],contact_matrix[7,5],
contact_matrix[7,6],contact_matrix[7,7],contact_matrix[7,8],contact_matrix[7,9],contact_matrix[7,10],
contact_matrix[8,1],contact_matrix[8,2],contact_matrix[8,3],contact_matrix[8,4],contact_matrix[8,5],
contact_matrix[8,6],contact_matrix[8,7],contact_matrix[8,8],contact_matrix[8,9],contact_matrix[8,10],
contact_matrix[9,1],contact_matrix[9,2],contact_matrix[9,3],contact_matrix[9,4],contact_matrix[9,5],
contact_matrix[9,6],contact_matrix[9,7],contact_matrix[9,8],contact_matrix[9,9],contact_matrix[9,10],
contact_matrix[10,1],contact_matrix[10,2],contact_matrix[10,3],contact_matrix[10,4],contact_matrix[10,5],
contact_matrix[10,6],contact_matrix[10,7],contact_matrix[10,8],contact_matrix[10,9],contact_matrix[10,10]};
int x_i[10] = {N_data[1],N_data[2],N_data[3],N_data[4],N_data[5],N_data[6],N_data[7],N_data[8],N_data[9],N_data[10]};
}
// ----------------------------------------------------------------------------------
parameters {
real<lower=0> phi_inv_0_9;
real<lower=0> phi_inv_10_19;
real<lower=0> phi_inv_20_29;
real<lower=0> phi_inv_30_39;
real<lower=0> phi_inv_40_49;
real<lower=0> phi_inv_50_59;
real<lower=0> phi_inv_60_69;
real<lower=0> phi_inv_70_79;
real<lower=0> phi_inv_80_89;
real<lower=0> phi_inv_over90;
real<lower=0> beta;
real<lower=0> gamma_0_9;
real<lower=0> gamma_10_19;
real<lower=0> gamma_20_29;
real<lower=0> gamma_30_39;
real<lower=0> gamma_40_49;
real<lower=0> gamma_50_59;
real<lower=0> gamma_60_69;
real<lower=0> gamma_70_79;
real<lower=0> gamma_80_89;
real<lower=0> gamma_over90;
}
// ----------------------------------------------------------------------------
transformed parameters{
real y_ode[n_days, 30]; // 30 because 10 of S + 10 of I + 10 of R (40_49 and 50_59)
real phi_0_9 = 1. / phi_inv_0_9;
real phi_10_19 = 1. / phi_inv_10_19;
real phi_20_29 = 1. / phi_inv_20_29;
real phi_30_39 = 1. / phi_inv_30_39;
real phi_40_49 = 1. / phi_inv_40_49;
real phi_50_59 = 1. / phi_inv_50_59;
real phi_60_69 = 1. / phi_inv_60_69;
real phi_70_79 = 1. / phi_inv_70_79;
real phi_80_89 = 1. / phi_inv_80_89;
real phi_over90 = 1. / phi_inv_over90;
{
real theta[11];
theta[1] = beta;
theta[2] = gamma_0_9;
theta[3] = gamma_10_19;
theta[4] = gamma_20_29;
theta[5] = gamma_30_39;
theta[6] = gamma_40_49;
theta[7] = gamma_50_59;
theta[8] = gamma_60_69;
theta[9] = gamma_70_79;
theta[10] = gamma_80_89;
theta[11] = gamma_over90;
y_ode = integrate_ode_rk45(sir, y0, initial_time, ts, theta, x_r, x_i);
}
}
// ----------------------------------------------------------------------------
model {
//priors
beta ~ normal(4, 0.1); //truncated at 0
gamma_0_9 ~ normal(0.5, 0.1); //truncated at 0
gamma_10_19 ~ normal(0.5, 0.1); //truncated at 0
gamma_20_29 ~ normal(0.5, 0.1); //truncated at 0
gamma_30_39 ~ normal(0.5, 0.1); //truncated at 0
gamma_40_49 ~ normal(0.5, 0.1); //truncated at 0
gamma_50_59 ~ normal(0.5, 0.1); //truncated at 0
gamma_60_69 ~ normal(0.5, 0.1); //truncated at 0
gamma_70_79 ~ normal(0.5, 0.1); //truncated at 0
gamma_80_89 ~ normal(0.5, 0.1); //truncated at 0
gamma_over90 ~ normal(0.5, 0.1); //truncated at 0
phi_inv_0_9 ~ exponential(5);
phi_inv_10_19 ~ exponential(5);
phi_inv_20_29 ~ exponential(5);
phi_inv_30_39 ~ exponential(5);
phi_inv_40_49 ~ exponential(5);
phi_inv_50_59 ~ exponential(5);
phi_inv_60_69 ~ exponential(5);
phi_inv_70_79 ~ exponential(5);
phi_inv_80_89 ~ exponential(5);
phi_inv_over90 ~ exponential(5);
//sampling distribution
//col(matrix x, int n) - The n-th column of matrix x. Here the number of infected people
// if (compute_likelihood == 1){
// I_0_9_data ~ neg_binomial_2(col(to_matrix(y_ode), 2) ,phi_0_9 );
// I_10_19_data ~ neg_binomial_2(col(to_matrix(y_ode), 5) ,phi_10_19);
// I_20_29_data ~ neg_binomial_2(col(to_matrix(y_ode), 8) ,phi_10_19);
// I_30_39_data ~ neg_binomial_2(col(to_matrix(y_ode), 11),phi_10_19);
// I_40_49_data ~ neg_binomial_2(col(to_matrix(y_ode), 14),phi_10_19);
// I_50_59_data ~ neg_binomial_2(col(to_matrix(y_ode), 17),phi_10_19);
// I_60_69_data ~ neg_binomial_2(col(to_matrix(y_ode), 20),phi_10_19);
// I_70_79_data ~ neg_binomial_2(col(to_matrix(y_ode), 23),phi_10_19);
// I_80_89_data ~ neg_binomial_2(col(to_matrix(y_ode), 26),phi_10_19);
// I_over90_data ~ neg_binomial_2(col(to_matrix(y_ode), 29),phi_10_19);
// }
}
// -------------------------------------------
generated quantities {
real recovery_time_0_9 = 1 / gamma_0_9;
real recovery_time_10_19= 1 / gamma_10_19;
real recovery_time_20_29 = 1 / gamma_20_29;
real recovery_time_30_39= 1 / gamma_30_39;
real recovery_time_40_49 = 1 / gamma_40_49;
real recovery_time_50_59= 1 / gamma_50_59;
real recovery_time_60_69 = 1 / gamma_60_69;
real recovery_time_70_79= 1 / gamma_70_79;
real recovery_time_80_89 = 1 / gamma_80_89;
real recovery_time_9over = 1 / gamma_over90;
real pred_I_0_9[n_days];
real pred_I_10_19[n_days];
real pred_I_20_29[n_days];
real pred_I_30_39[n_days];
real pred_I_40_49[n_days];
real pred_I_50_59[n_days];
real pred_I_60_69[n_days];
real pred_I_70_79[n_days];
real pred_I_80_89[n_days];
real pred_I_over90[n_days];
pred_I_0_9 = neg_binomial_2_rng(col(to_matrix(y_ode), 2) + 1e-5, phi_0_9);
pred_I_10_19 = neg_binomial_2_rng(col(to_matrix(y_ode), 5) + 1e-5, phi_10_19);
pred_I_20_29 = neg_binomial_2_rng(col(to_matrix(y_ode), 8) + 1e-5, phi_20_29);
pred_I_30_39 = neg_binomial_2_rng(col(to_matrix(y_ode), 11) + 1e-5, phi_30_39);
pred_I_40_49 = neg_binomial_2_rng(col(to_matrix(y_ode), 14) + 1e-5, phi_40_49);
pred_I_50_59 = neg_binomial_2_rng(col(to_matrix(y_ode), 17) + 1e-5, phi_50_59);
pred_I_60_69 = neg_binomial_2_rng(col(to_matrix(y_ode), 20) + 1e-5, phi_60_69);
pred_I_70_79 = neg_binomial_2_rng(col(to_matrix(y_ode), 23) + 1e-5, phi_70_79);
pred_I_80_89 = neg_binomial_2_rng(col(to_matrix(y_ode), 26) + 1e-5, phi_80_89);
pred_I_over90 = neg_binomial_2_rng(col(to_matrix(y_ode), 29) + 1e-5, phi_over90);
}