Hello again!

@wds15, @yizhang and @bbbales2, per @maxbiostat 's suggestion, I am seeking your advice for my question. I’ve given it my best try, and what follows is as far as I’ve gotten on this problem:

```
functions {
// This largely follows the deSolve package, but also includes the x_r and x_i variables.
// These variables are used in the background.
real[] deltE(real time,
real[] peakE,
real[] wrE,
real[] wfE) {
for (i in 1:size(t)){
return t<peakE ? exp((-(t-peakE)^2)/wrE) : exp((-(t-peakE))/wfE);
}
}
real[] deltG(real time,
real[] peakG,
real[] wrG,
real[] wfG) {
for (i in 1:size(t)){
return t<peakG ? exp((-(t-peakG)^2)/wrG) : exp((-(t-peakG))/wfG);
}
}
real[] deltM(real time,
real[] peakM,
real[] wrM,
real[] wfM) {
for (i in 1:size(t)){
return t<peakM ? exp((-(t-peakM)^2)/wrM) : exp((-(t-peakM))/wfM);
}
}
real[] TIV(real time,
real[] y, // system components {target cell, infected cell, viral load}
real[] params, // parameters
real[] x_r,
int[] x_i) {
real T = y[1];
real I = y[2];
real V = y[3];
real N = x_i[1];
real d = params[1];
real delta = params[2];
real c = params[3];
real kt.Aa = params[4];
real p = params[5];
real ktauE = params[6];
real ktauG = params[7];
real ktauM = params[8];
for (i in 1:t) {
dT_dt = (params[1] * 1e+12) - (params[1] * y[1]) - (params[4] * y[1] * y[3]);
dI_dt = (params[4] * y[1] * y[3]) - ((params[2] + params[6] * deltE[i]) * y[2]);
dV_dt = (params[5] * y[2]) - ((params[3] + params[7] * deltG[i] + params[8] * deltM[i]) * y[3]);
return {dT_dt, dI_dt, dV_dt};
}
}
}
data {
int<lower = 1> n_obs; // Number of obs
real y0[3];
real t0;
int<lower = 1> n_params; // Number of model parameters
int<lower = 1> n_difeq; // Number of differential equations in the system
int<lower = 1> n_fake; // This is to generate "predicted"/"unsampled" data
real ts[n_obs]; // Time points that were sampled
int N;
real fake_ts[n_fake]; // Time points for "predicted"/"unsampled" data
}
transformed data {
real x_r[0];
int x_i[1] = {N};
}
parameters {
real<lower = 0> params[n_params]; // Model parameters
real<lower = 0> V0; // Initial viral load
}
transformed parameters{
real y_hat[n_obs, n_difeq]; // Output from the ODE solver
real y0[n_difeq]; // Initial conditions for both S and I
{
real params[8];
real params[1] = d;
real params[2] = delta;
real params[3] = c;
real params[4] = kt.Aa;
real params[5] = p;
real params[6] = ktauE;
real params[7] = ktauG;
real params[8] = ktauM;
y_hat = integrate_ode_rk45(TIV, y0, t0, ts, params, x_r, x_i);
}
}
model {
//priors
d ~ normal(0, 0.5); //truncated at 0
delta ~ normal(0, 1);
c ~ normal(0, 1);
kt.Aa ~ normal(0, 1);
p ~ normal(0, 1);
ktauE ~ normal(0, 1);
ktauG ~ normal(0, 1);
ktauM ~ normal(0, 1);
//sampling distribution
//col(matrix x, int n) - The n-th column of matrix x. Here the number of infected people
cases ~ neg_binomial_2(col(to_matrix(y), 2), phi);
}
generated quantities {
// Generate predicted data over the whole time series:
real fake_I[n_fake, n_difeq];
fake_I = integrate_ode_rk45(SI, y0, t0, fake_ts, params, x_r, x_i);
}
```

Again, many thanks for any help offered!

EDIT: @maxbiostat edited this post for syntax highlighting.