# Explaining how Rstan estimates parameter using ODE and data.!

**URL:** <https://discourse.mc-stan.org/t/explaining-how-rstan-estimates-parameter-using-ode-and-data/21478>\
**Category:** General\
**Created:** [March 24, 2021, 11:50pm UTC](https://discourse.mc-stan.org/t/explaining-how-rstan-estimates-parameter-using-ode-and-data/21478 "2021-03-24T23:50:26Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Boyh](https://avatars.discourse-cdn.com/v4/letter/b/a9adbd/32.png) [@Boyh](https://discourse.mc-stan.org/u/Boyh)\
**Post date:** [March 24, 2021, 11:50pm UTC](https://discourse.mc-stan.org/t/explaining-how-rstan-estimates-parameter-using-ode-and-data/21478/1 "2021-03-24T23:50:27Z")

</div>

Please can anyone be kind enough to explain the process that Rstan uses to estimate parameters given a differential equation and collected data. A link to a page or journal article will be appreciated. Thanks

---

<div class="post-metadata">

**Author:** ![wds15](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/wds15/32/908_2.png) [@wds15](https://discourse.mc-stan.org/u/wds15)\
**Post date:** [March 25, 2021, 8:09am UTC](https://discourse.mc-stan.org/t/explaining-how-rstan-estimates-parameter-using-ode-and-data/21478/2 "2021-03-25T08:09:50Z")

</div>

Like here:

- [Bayesian Model of Planetary Motion: exploring ideas for a modeling workflow when dealing with ordinary differential equations and multimodality | planetary\_motion.utf8](https://mc-stan.org/users/documentation/case-studies/planetary_motion/planetary_motion.html)

- [Upgrading to the new ODE interface](https://mc-stan.org/users/documentation/case-studies/convert_odes.html)

- [Bayesian workflow for disease transmission modeling in Stan](https://mc-stan.org/users/documentation/case-studies/boarding_school_case_study.html)

?

---

<div class="post-metadata">

**Author:** ![srossell](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/srossell/32/10567_2.png) [@srossell](https://discourse.mc-stan.org/u/srossell)\
**Post date:** [March 25, 2021, 8:25am UTC](https://discourse.mc-stan.org/t/explaining-how-rstan-estimates-parameter-using-ode-and-data/21478/3 "2021-03-25T08:25:42Z")

</div>

Hi,

I wrote a blogpost with a very simple example.  
[https://towardsdatascience.com/bayesian-inference-and-differential-equations-cb646b50f148](https://towardsdatascience.com/bayesian-inference-and-differential-equations-cb646b50f148)

I hope it will be useful for you

---

<div class="post-metadata">

**Author:** ![Boyh](https://avatars.discourse-cdn.com/v4/letter/b/a9adbd/32.png) [@Boyh](https://discourse.mc-stan.org/u/Boyh)\
**Post date:** [March 27, 2021, 12:57pm UTC](https://discourse.mc-stan.org/t/explaining-how-rstan-estimates-parameter-using-ode-and-data/21478/4 "2021-03-27T12:57:48Z")

</div>

Thank you @srossell , thank you @wds15, I have been working with [Bayesian workflow for disease transmission modeling in Stan](https://mc-stan.org/users/documentation/case-studies/boarding_school_case_study.html) to estimate the SIR model parameter values on Covid of a country using the same priors as used in the work , but I keep getting the messages belows

"1: The largest R-hat is 1.73, indicating chains have not mixed.  
Running the chains for more iterations may help. See  
[http://mc-stan.org/misc/warnings.html#r-hat](http://mc-stan.org/misc/warnings.html#r-hat)  
2: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.  
Running the chains for more iterations may help. See  
[http://mc-stan.org/misc/warnings.html#bulk-ess](http://mc-stan.org/misc/warnings.html#bulk-ess)  
3: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.  
Running the chains for more iterations may help. See  
[http://mc-stan.org/misc/warnings.html#tail-ess](http://mc-stan.org/misc/warnings.html#tail-ess)  
"  
I ran the chain up to 7000 iterations and visited the suggested links, yet the chains are still not mixing. I am a maths student however and I am trying to master Bayern and Rstan but is really not working. Thanks

---

<div class="post-metadata">

**Author:** ![srossell](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/srossell/32/10567_2.png) [@srossell](https://discourse.mc-stan.org/u/srossell)\
**Post date:** [March 28, 2021, 12:35pm UTC](https://discourse.mc-stan.org/t/explaining-how-rstan-estimates-parameter-using-ode-and-data/21478/5 "2021-03-28T12:35:08Z")

</div>

Hi Boyh,

Your question is too open. If you are using Grinsztajn et al’s Stan program, it may be best to approach one of the authors. If you’ve edited their code, then make incremental changes to Grinsztajn et al’s code. Hopefully, in that way you’ll identify you issue.

Good luck

---

<div class="post-metadata">

**Author:** ![Funko\_Unko](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/funko_unko/32/10786_2.png) [@Funko\_Unko](https://discourse.mc-stan.org/u/Funko_Unko)\
**Post date:** [March 28, 2021, 9:39pm UTC](https://discourse.mc-stan.org/t/explaining-how-rstan-estimates-parameter-using-ode-and-data/21478/6 "2021-03-28T21:39:19Z")

</div>

Hmmmm, that paper / model was written at the beginning of the pandemic, wasn’t it?

Maybe it works better with early data, and has to be extended to work up to current events?

---

<div class="post-metadata">

**Author:** ![Boyh](https://avatars.discourse-cdn.com/v4/letter/b/a9adbd/32.png) [@Boyh](https://discourse.mc-stan.org/u/Boyh)\
**Post date:** [March 29, 2021, 8:42am UTC](https://discourse.mc-stan.org/t/explaining-how-rstan-estimates-parameter-using-ode-and-data/21478/7 "2021-03-29T08:42:39Z")

</div>

yes @Funko_Unko It is and I also using it for the data at the beginning of the pandemic., but I keep getting the followings below

1: The largest R-hat is 1.73, indicating chains have not mixed.  
Running the chains for more iterations may help. See  
[http://mc-stan.org/misc/warnings.html#r-hat](http://mc-stan.org/misc/warnings.html#r-hat)  
2: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.  
Running the chains for more iterations may help. See  
[http://mc-stan.org/misc/warnings.html#bulk-ess](http://mc-stan.org/misc/warnings.html#bulk-ess)  
3: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.  
Running the chains for more iterations may help. See  
[http://mc-stan.org/misc/warnings.html#tail-ess](http://mc-stan.org/misc/warnings.html#tail-ess)

I used the prior model used, adjusted it and even increased the iteration number to 7000 (just for 216 observable data), but the chains ain’t still mixing well. R hat used for determining the accuracy of the whole process is still greater than 1.  
I do not know what else to do. Thanks

---

<div class="post-metadata">

**Author:** ![maxbiostat](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/maxbiostat/32/1326_2.png) [@maxbiostat](https://discourse.mc-stan.org/u/maxbiostat)\
**Post date:** [March 29, 2021, 1:21pm UTC](https://discourse.mc-stan.org/t/explaining-how-rstan-estimates-parameter-using-ode-and-data/21478/8 "2021-03-29T13:21:29Z")

</div>

Let’s back up a bit. Can you post your model and perhaps some data so we can have a look? The resources @wds15 posted are also worth checking out.

---

<div class="post-metadata">

**Author:** ![Boyh](https://avatars.discourse-cdn.com/v4/letter/b/a9adbd/32.png) [@Boyh](https://discourse.mc-stan.org/u/Boyh)\
**Post date:** [March 29, 2021, 3:45pm UTC](https://discourse.mc-stan.org/t/explaining-how-rstan-estimates-parameter-using-ode-and-data/21478/9 "2021-03-29T15:45:30Z")

</div>

Attached below is the data and codes. Thank you  
[covid\_data.csv](https://discourse.mc-stan.org/uploads/short-url/4WAQ1bE62IjKfhvclS3xbNXKAgy.csv) (3.3 KB)

```
//saved as seir.stan
functions {
  real[] sir(real t, real[] y, real[] theta, 
             real[] x_r, int[] x_i) {

      real S = y[1];
      real I = y[2];
      real R = y[3];
      int N = x_i[1];
      
      real beta = theta[1];
      real gamma = theta[2];
      
      real dS_dt = -beta * I * S / N;
      real dI_dt = beta * I * S / N - gamma * I;
      real dR_dt = gamma * I;
      
      return {dS_dt, dI_dt, dR_dt};
  }
}
data {
  int<lower=1> n_days;
  real t0;
  real y0[3];
  real ts[n_days];
  int N;
  int cases[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> phi_inv;
  real<lower=0, upper=1> p_reported; // proportion of infected (symptomatic) people reported
}
transformed parameters{
  real y[n_days, 3];
  real incidence[n_days - 1];
  real phi = 1. / phi_inv;
  // initial compartement values
  real theta[2] = {beta, gamma};
  y = integrate_ode_rk45(sir, y0, t0, ts, theta, x_r, x_i);
  for (i in 1:n_days-1){
    incidence[i] = (y[i, 1] - y[i+1, 1]) * p_reported;
  }
}
model {
  //priors
  beta ~ normal(2, 1);
  gamma ~ normal(0.9, 1);
  phi_inv ~ exponential(5);
  p_reported ~ beta(1, 2);
  //sampling distribution
  //col(matrix x, int n) - The n-th column of matrix x. Here the number of infected people 
  cases[1:(n_days-1)] ~ neg_binomial_2(incidence, phi);
}
generated quantities {
  real R0 = beta / gamma;
  real recovery_time = 1 / gamma;
  real pred_cases[n_days-1];
  pred_cases = neg_binomial_2_rng(incidence, phi);
}

```

```
setwd("~/R")

nga <- read.csv("covid_data.csv", header=TRUE, sep=",")

# Nigeria population
N <- 211400708;

# time series of cases
cases <- nga$Number_of_Reported_Cases # Number of cases

# times
n_days <- length(cases) 
t <- seq(0, n_days, by = 1)
t0 = 0 
t <- t[-1]

#initial conditions
i0 <- 5
s0 <- N - i0
r0 <- 0
y0 = c(S = s0, I = i0, R = r0)

# data for Stan
data_sir <- list(n_days = n_days, y0 = y0, t0 = t0, ts = t, N = N, cases = cases)

# number of MCMC steps
niter <- 2000

model <- stan_model("seir.stan")

#pass the data to stan and run the model
fit_seir <- sampling(model,
                           data = data_sir,
                           iter = niter,
                           chains = 4, 
                           seed = 0,
                           control = list(max_treedepth = 15, adapt_delta=0.99))

traceplot(fit_seir, pars = c('beta', 'gamma', "R0", "recovery_time"))
pars=c('beta', 'gamma', "R0", "recovery_time")
print(fit_seir, pars = pars)
stan_dens(fit_seir, pars = pars, separate_chains = TRUE)

```

---

<div class="post-metadata">

**Author:** ![maxbiostat](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/maxbiostat/32/1326_2.png) [@maxbiostat](https://discourse.mc-stan.org/u/maxbiostat)\
**Post date:** [March 29, 2021, 6:52pm UTC](https://discourse.mc-stan.org/t/explaining-how-rstan-estimates-parameter-using-ode-and-data/21478/10 "2021-03-29T18:52:19Z")

</div>

Thanks. I’m a bit swapped for the coming days, so will tag @charlesm93 just in case he has some free time.

---

<div class="post-metadata">

**Author:** ![Boyh](https://avatars.discourse-cdn.com/v4/letter/b/a9adbd/32.png) [@Boyh](https://discourse.mc-stan.org/u/Boyh)\
**Post date:** [March 31, 2021, 9:08pm UTC](https://discourse.mc-stan.org/t/explaining-how-rstan-estimates-parameter-using-ode-and-data/21478/11 "2021-03-31T21:08:12Z")

</div>

@charlesm93 , I know you must be really busy, please can you spend few mins and look into my codes, it really matters to me. I have done everything I can. Thank you.

---

<div class="post-metadata">

**Author:** ![charlesm93](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/charlesm93/32/4634_2.png) [@charlesm93](https://discourse.mc-stan.org/u/charlesm93)\
**Post date:** [April 1, 2021, 12:50am UTC](https://discourse.mc-stan.org/t/explaining-how-rstan-estimates-parameter-using-ode-and-data/21478/12 "2021-04-01T00:50:41Z")

</div>

Hi,  
Can you print the trace plots, including the warmup phase, and the summary of the fit? Please include `lp__` in the parameters, that is

```
pars = c('beta', 'gamma', 'R0', 'recovery_time', 'lp__')

```

Once I see these, I’ll make further comments.

---

<div class="post-metadata">

**Author:** ![maxbiostat](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/maxbiostat/32/1326_2.png) [@maxbiostat](https://discourse.mc-stan.org/u/maxbiostat)\
**Post date:** [April 1, 2021, 12:54am UTC](https://discourse.mc-stan.org/t/explaining-how-rstan-estimates-parameter-using-ode-and-data/21478/13 "2021-04-01T00:54:45Z")

</div>

I can.

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/7/71cc121fb1a3aeeb168a5ce31b5b601afe4757c9.jpeg)

I used cmdstanr running cmdstan-2.26.1 though.

EDIT: And I further excluded a point (cases[107]) which was suspiciously 0 in the middle of the epidemic.

---

<div class="post-metadata">

**Author:** ![charlesm93](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/charlesm93/32/4634_2.png) [@charlesm93](https://discourse.mc-stan.org/u/charlesm93)\
**Post date:** [April 1, 2021, 12:57am UTC](https://discourse.mc-stan.org/t/explaining-how-rstan-estimates-parameter-using-ode-and-data/21478/14 "2021-04-01T00:57:31Z")

</div>

Thanks, can you also include the warmup phase?

---

<div class="post-metadata">

**Author:** ![maxbiostat](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/maxbiostat/32/1326_2.png) [@maxbiostat](https://discourse.mc-stan.org/u/maxbiostat)\
**Post date:** [April 1, 2021, 1:03am UTC](https://discourse.mc-stan.org/t/explaining-how-rstan-estimates-parameter-using-ode-and-data/21478/15 "2021-04-01T01:03:19Z")

</div>

Sorry, you did ask for that and I forgot. Here’s what I get when using `inc_warmup=TRUE`.

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/1/1667eff80af3239549d8544227a64bc1eaf2a797.jpeg)

Looks weird. What gives?

---

<div class="post-metadata">

**Author:** ![charlesm93](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/charlesm93/32/4634_2.png) [@charlesm93](https://discourse.mc-stan.org/u/charlesm93)\
**Post date:** [April 1, 2021, 1:06am UTC](https://discourse.mc-stan.org/t/explaining-how-rstan-estimates-parameter-using-ode-and-data/21478/16 "2021-04-01T01:06:38Z")

</div>

I think you didn’t use the save\_warmup option if you fitted the model with cmdstanr. So it’s returning the same iterations.

---

<div class="post-metadata">

**Author:** ![maxbiostat](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/maxbiostat/32/1326_2.png) [@maxbiostat](https://discourse.mc-stan.org/u/maxbiostat)\
**Post date:** [April 1, 2021, 1:22am UTC](https://discourse.mc-stan.org/t/explaining-how-rstan-estimates-parameter-using-ode-and-data/21478/17 "2021-04-01T01:22:52Z")

</div>

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/1/16dcec6f7c4c8efccd627fa2d6fe7a7abb7d8111.jpeg)  
Used a few more iterations this time.

Seems like a case of the unidentifiability-itis.  
EDIT: Wrong diagnosis. Could be a case of multimodality-itis, but we’re waiting on the lab results [OK, I’m done with the disease analogy now].

---

<div class="post-metadata">

**Author:** ![Funko\_Unko](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/funko_unko/32/10786_2.png) [@Funko\_Unko](https://discourse.mc-stan.org/u/Funko_Unko)\
**Post date:** [April 1, 2021, 8:16am UTC](https://discourse.mc-stan.org/t/explaining-how-rstan-estimates-parameter-using-ode-and-data/21478/18 "2021-04-01T08:16:37Z")

</div>

For unidentifiability we would like the posterior predictive plots to look similar though, don’t we? Seems unlikely that they do.

---

<div class="post-metadata">

**Author:** ![maxbiostat](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/maxbiostat/32/1326_2.png) [@maxbiostat](https://discourse.mc-stan.org/u/maxbiostat)\
**Post date:** [April 1, 2021, 11:06am UTC](https://discourse.mc-stan.org/t/explaining-how-rstan-estimates-parameter-using-ode-and-data/21478/19 "2021-04-01T11:06:56Z")

</div>

They lead to the same R\_0, which ultimately controls the dynamics. Here are traceplots of selected points of the posterior predictive:

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/9/9ed32ca6c8ac072ac9edd1c0cd50bb8c8ea667a7.jpeg)

---

<div class="post-metadata">

**Author:** ![Funko\_Unko](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/funko_unko/32/10786_2.png) [@Funko\_Unko](https://discourse.mc-stan.org/u/Funko_Unko)\
**Post date:** [April 1, 2021, 11:19am UTC](https://discourse.mc-stan.org/t/explaining-how-rstan-estimates-parameter-using-ode-and-data/21478/20 "2021-04-01T11:19:32Z")

</div>

I can see exactly nothing on the latter two plots.

It does look though like that one chain consistently predicts much fewer cases for day 1 and day 20, where it should be in the order of zero.

[Next page](https://discourse.mc-stan.org/t/explaining-how-rstan-estimates-parameter-using-ode-and-data/21478.md?page=2)