# How to include a changing coefficient in ODEs in rstan

Hi all, I’m new to stan. I now have trouble in including a temperature-denpenent coefficient in stan. Below is my code. In this ODE, the coefficent “b” can be calculated based on a given temperature. However, I only know how to build a stan model by using a fixed value, but I have no idea on how to include a varying coefficent into stan.

``````functions {
real[] seir(real t, real[] y, real[] theta,
real[] x_r, int[] x_i) {

real N = x_i[1];

real beta = theta[1];
real gamma = theta[2];
real a = theta[3];
real i0 = theta[4];
real e0 = theta[5];
real b;

real init[4] = {N - i0 - e0, e0, i0, 0};
real S = y[1] + init[1];
real E = y[2] + init[2];
real I = y[3] + init[3];
real R = y[4] + init[4];

real dS_dt = -b * beta * I * S / N;
real dE_dt =  b * beta * I * S / N - a * E;
real dI_dt = a * E - gamma * I;
real dR_dt =  gamma * I;

return {dS_dt, dE_dt, dI_dt, dR_dt};
}
}
data {
int<lower=1> n_days;
real t0;
real y0[4];
real ts[n_days];
int N;
int cases[n_days];
real b[n_days];
}
transformed data {
real x_r[0];
int x_i[1] = { N };
}
parameters {
real<lower=0> gamma;
real<lower=0> beta;
real<lower=0> a;
real<lower=0> phi_inv;
real<lower=0> i0;
real<lower=0> e0;
}
transformed parameters{
real y[n_days, 4];
real phi = 1. / phi_inv;
real theta[5] = {beta, gamma, a, i0, e0};
y = integrate_ode_rk45(seir, rep_array(0.0, 4), t0, ts, theta, x_r, x_i);
}
model {
beta ~ normal(0,1);
gamma ~ normal(0,1);
a ~ normal(0,1);
phi_inv ~ exponential(5);
i0 ~ normal(0, 10);
e0 ~ normal(0, 10);

cases[1:(n_days)] ~ neg_binomial_2(y[,3], phi);
}
generated quantities {
real M = beta / gamma;
}

``````

My data looks like:
cases<-c(2,5,3,8,12,10,14, 20, 23…) (n=104)
b<-c(0.798, 0.752, 0.813, 0.806,0.774, 0.785,…) (n=104)

If I run the above code, I’ll get the error: “Exception: neg_binomial_2_lpmf: Location parameter[1] is nan, but must be > 0!” I know it must be something wrong with “b”, but I dont know how to correct it. Thank you so much if anyone can help me.

Hi Ruby,

Since b is not a parameter, you could pass b in x_r in the transformed data block. Then, inside your ode function, have a for statement which loops through time, combined with an if statement to assign to b the respective entry from x_r if the solver is currently at time point “t”.

Hope I got your question right.