# Turing package(s) + parameter estimation for DiffEq models?

**URL:** <https://discourse.julialang.org/t/turing-package-s-parameter-estimation-for-diffeq-models/44082>\
**Category:** New to Julia\
**Created:** [August 1, 2020, 11:10am UTC](https://discourse.julialang.org/t/turing-package-s-parameter-estimation-for-diffeq-models/44082 "2020-08-01T11:10:53Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [August 1, 2020, 11:10am UTC](https://discourse.julialang.org/t/turing-package-s-parameter-estimation-for-diffeq-models/44082/1 "2020-08-01T11:10:53Z")

</div>

Are there any illustrative notebooks around on how to do proper model fitting of dynamic models with parameter statistics in Turing tools? Can it be done?

To be more concrete… take the SIR epidemiology model below with boarding school data for measles in North England in 1978 – data and model parameters are taken from Maia Martcheva’s Springer book _An Introduction to Mathematical Epidemiology_ from 2015.

```julia
# Packages used
using Plots; pyplot()
using LaTeXStrings
using DifferentialEquations
#
N = 763 # Number of pupils in boarding school
tm = 3:14 # days of infection data
# Infection data
Id = [25, 75, 227, 296, 258, 236, 192, 126, 71, 28, 11, 7]

```

SIR model function – S, I, R, N are _number of persons_ (i.e., not per capita) – “my” parameters are \tau\_\mathrm{i} for infection time (equal to \tau\_\mathrm{i} = 1/(N\cdot \beta) in Martcheva’s book) and \tau\_\mathrm{r} for recovery time (equal to \tau\_\mathrm{r} = 1/\alpha in Martcheva’s book, or 1/\mu in some other books). See p. 126-129 in Martcheva’s book for data and parameters. The function is:

```julia
# Function for SIR model
function sir!(dx,x,p,t)
    S,I,R = x
    tau_i,tau_r,N = p
    #
    dS = -I*S/tau_i/N
    dI = I*S/tau_i/N - I/tau_r
    dR = I/tau_r
    dx .= [dS,dI,dR]
end
;

```

Problem set-up:

```julia
I3 = Float64(Id[1])
S3 = N-I3
R3 = 0.
x3 = [S3,I3,R3]
#
tau_i = 0.553
tau_r = 2.15
RR0 = round(tau_r/tau_i,digits=1) # The "R0" parameter
p = [tau_i, tau_r,N]
#
tspan=(3.,14.)
#
prob = ODEProblem(sir!,x3,tspan,p)
sol = solve(prob)
;

```

Parameters for model are suggested after some model fitting in Martcheva’s book. The results are as follows:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/4/c/4c15630330b34f7f8dffdd188ad6cf1fac04dc4b.png)

_Question_: Suppose I want to _redo_ the model-fit, and also find some description of model uncertainty for the model parameters… Can this be done using Turing? Any notebook examples?

In this particular case, the initial states are probably well known. But in general, one would also want to estimate those.

Any suggestions and help are appreciated.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 1, 2020, 12:38pm UTC](https://discourse.julialang.org/t/turing-package-s-parameter-estimation-for-diffeq-models/44082/2 "2020-08-01T12:38:12Z")

</div>

[https://turing.ml/dev/tutorials/10-bayesiandiffeq/](https://turing.ml/dev/tutorials/10-bayesiandiffeq/)

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [August 5, 2020, 9:38am UTC](https://discourse.julialang.org/t/turing-package-s-parameter-estimation-for-diffeq-models/44082/3 "2020-08-05T09:38:25Z")

</div>

Thanks for link – I’m trying to work my way through it. I have now installed Julia v. 1.5.0, and just updated all packages.

First problem:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/3/5/35fb20c70935c61667426f0e8b14ae85b201be31.png)

There is a problem with package `Turing`, or perhaps with package `MCMCChains` – `Turing` doesn’t find `MCMCChains.AbstractChains`.

Also, from the linked document: it is not 100% clear which packages I use where. Do I need to import `Turing` in order to run `DiffEqBayes`, etc.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [August 5, 2020, 9:53am UTC](https://discourse.julialang.org/t/turing-package-s-parameter-estimation-for-diffeq-models/44082/4 "2020-08-05T09:53:21Z")

</div>

> [@ChrisRackauckas](#):
>
> [https://turing.ml/dev/tutorials/10-bayesiandiffeq/](https://turing.ml/dev/tutorials/10-bayesiandiffeq/)

In addition to the problem of importing `Turing`, using `DiffEqBayes`, etc. due to error messages:

1. In section “Fitting Lotka-Volterra with DiffEqBayes”, suddenly one parameter shows up which has not been declared: \sigma. It is not clear from the context what parameter \sigma represents.
2. If I “cheat” and look at the _next_ section, “Direct Handling of Bayesian Estimation with Turing”, in that case, parameter \sigma is defined explicitly. And it is proposed as a prior of the unknown _variance_ (or _covariance_ if you will) of the normally distributed noise in the data (the second parameter in `MvNormal` is variance/covariance, while strangely, the second parameter in `Normal` is standard deviation).
3. From Wikipedia, it is common in Bayes estimation to assume an `InverseGamma` distribution for a parameter which represents _variance_. (It is not clear to me why parameters (2,3) are chosen.)

Anyway, it makes sense that \sigma is the variance of the (normally distributed) observation noise.

**Question** : is my assumption correct? In other words, when I use `DiffEqBayes`, parameter \sigma is _automatically_ added?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 5, 2020, 9:56am UTC](https://discourse.julialang.org/t/turing-package-s-parameter-estimation-for-diffeq-models/44082/5 "2020-08-05T09:56:28Z")

</div>

> [@BLI](#):
>
> Anyway, it makes sense that σ\sigma is the variance of the (normally distributed) observation noise.
> 
> **Question** : is my assumption correct? In other words, when I use `DiffEqBayes` , parameter σ\sigma is _automatically_ added?

Yes, @Vaibhavdixit02 we should make sure this is a bit more clearly explained in the tutorial. It is measurement noise on the outputs, it’s common to use `InverseGamma`, those parameters are just “a good prior” (tends to work, just fairly positive with non-zero variance and moves easily), but it’s also just a choice.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [August 5, 2020, 10:01am UTC](https://discourse.julialang.org/t/turing-package-s-parameter-estimation-for-diffeq-models/44082/6 "2020-08-05T10:01:48Z")

</div>

> [@ChrisRackauckas](#):
>
> [https://turing.ml/dev/tutorials/10-bayesiandiffeq/](https://turing.ml/dev/tutorials/10-bayesiandiffeq/)

In section “Fitting Lotka-Volterra with DiffEqBayes”, plotting of the `bayesian_result_turing` leads to a 5\times 2 matrix of plots.

1. In the first plot column, there is one plot for each parameter and how it varies with “iteration” (which probably refers to `num_samples` in the specification of the `turing_inference` problem).  
Does “iteration” here refer to the iteration of the estimate towards a minimum (as in LS parameter fitting)? Or does it (more likely…) refer to the number of parameter estimates, much like the number of Monte Carlo estimates in Bootstrapping?

2. I _assume_ that the plots in column 2 – the parameter estimate distributions – are simply a “histograms” of the Sample value vs. iteration in column 1?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 5, 2020, 10:03am UTC](https://discourse.julialang.org/t/turing-package-s-parameter-estimation-for-diffeq-models/44082/7 "2020-08-05T10:03:01Z")

</div>

> [@BLI](#):
>
> Does “iteration” here refer to the iteration of the estimate towards a minimum (as in LS parameter fitting)? Or does it (more likely…) refer to the number of parameter estimates, much like the number of Monte Carlo estimates in Bootstrapping?

It’s the Monte Carlo chain, the MCMC choices. It’s the chain that over time should be well-mixed to give the distribution of the parameters.

> [@BLI](#):
>
> I _assume_ that the plots in column 2 – the parameter estimate distributions – are simply a “histograms” of the Sample value vs. iteration in column 1?

Exactly.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [August 5, 2020, 10:04am UTC](https://discourse.julialang.org/t/turing-package-s-parameter-estimation-for-diffeq-models/44082/8 "2020-08-05T10:04:47Z")

</div>

Excellent. If I now just could get `Turing` and the other packages working…

---

<div class="post-metadata">

**Author:** ![Vaibhavdixit02](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vaibhavdixit02/32/2916_2.png) [@Vaibhavdixit02](https://discourse.julialang.org/u/Vaibhavdixit02)\
**Post date:** [August 5, 2020, 10:14am UTC](https://discourse.julialang.org/t/turing-package-s-parameter-estimation-for-diffeq-models/44082/9 "2020-08-05T10:14:33Z")

</div>

@BLI we have been discussing about updating this tutorial on slack.

> @Vaibhavdixit02 we should make sure this is a bit more clearly explained in the tutorial.

Yes this should be updated.

> If I now just could get `Turing` and the other packages working…

I think you might be experiencing the same issue as [Slack](https://julialang.slack.com/archives/CCYDC34A0/p1596567897449200), you might want to try with Turing#master

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [August 5, 2020, 12:48pm UTC](https://discourse.julialang.org/t/turing-package-s-parameter-estimation-for-diffeq-models/44082/10 "2020-08-05T12:48:52Z")

</div>

Thanks! Now it works. (On my laptop.)

Question: I ran the “Fitting Lotka-Volterra with DiffEqBayes” example, leading to `bayesian_result_turing` which I then plotted.

**Question** : Is there a way to get out the data manually of this structure `bayesian_result_turing`? Specifically, can I extract the parameter distributions? And the Sample values? If I can get out the Sample values, I can plot them as scatter plots (Sample plot 1 vs. Sample plot 2), and get ideas of correlations in parameter estimates – I’d assume… Or produce 2D histograms/fit surfaces to the 2D histograms…

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [August 5, 2020, 2:26pm UTC](https://discourse.julialang.org/t/turing-package-s-parameter-estimation-for-diffeq-models/44082/11 "2020-08-05T14:26:02Z")

</div>

In Turing (and other languages), the \sim symbol appears…, e.g.  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/7/0/70ad3ef09a126be3b5bd5c46ad7d0cdf002b5b65.png)

However, the intro to Turing _fails_ to explain which symbol this is. I have tried `\sim` followed by TAB, and that gives a similar look, but `Turing` doesn’t understand this symbol.

So… which symbol should I use?

---

<div class="post-metadata">

**Author:** ![arnaudmgh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arnaudmgh/32/7134_2.png) [@arnaudmgh](https://discourse.julialang.org/u/arnaudmgh)\
**Post date:** [August 5, 2020, 3:07pm UTC](https://discourse.julialang.org/t/turing-package-s-parameter-estimation-for-diffeq-models/44082/12 "2020-08-05T15:07:31Z")

</div>

Hi, you can access the values with `bayesian_result_turing.value`.  
to get the accessible field names of a variable you can type: `fieldnames(typeof(bayesian_result_turing))`

I am in the process of updating the tutorial, with the help of Vaibhav, Chris and others, if you want to take a peak: [https://github.com/arnaudmgh/TuringTutorials/blob/master/10\_diffeq.ipynb](https://github.com/arnaudmgh/TuringTutorials/blob/master/10_diffeq.ipynb)

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [August 5, 2020, 3:09pm UTC](https://discourse.julialang.org/t/turing-package-s-parameter-estimation-for-diffeq-models/44082/13 "2020-08-05T15:09:02Z")

</div>

Do you also know which symbol I should use for assigning distributions… symbol similar to \sim, but different. Tutorials are not really readable if it is presumed that the reader knows how to type in such special symbols.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [August 5, 2020, 3:13pm UTC](https://discourse.julialang.org/t/turing-package-s-parameter-estimation-for-diffeq-models/44082/14 "2020-08-05T15:13:06Z")

</div>

It would also be nice if you include 1-2 examples where:

- only a subset of the states are observed (can `save_idx` be used?), and
- a nonlinear mapping of states are observed, e.g., y=g(x) where x\in \mathbb{R}^n, y\in \mathbb{R}^m, and normally (but not always?) n\ge m. (Would you have to use callbacks for this?)

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [August 5, 2020, 3:28pm UTC](https://discourse.julialang.org/t/turing-package-s-parameter-estimation-for-diffeq-models/44082/15 "2020-08-05T15:28:34Z")

</div>

Cool…:  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/2/a/2a5def3700661a27bc427bba40a4e120c6667fcc.png)

Is there a simple way to fit a surface to this histogram, and produce a contour plot?

Alternative, with `scatter` plot:  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/6/f/6fa209872d84c0784a5d857b8755fcb8013b04b5.png)

---

<div class="post-metadata">

**Author:** ![arnaudmgh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arnaudmgh/32/7134_2.png) [@arnaudmgh](https://discourse.julialang.org/u/arnaudmgh)\
**Post date:** [August 5, 2020, 3:37pm UTC](https://discourse.julialang.org/t/turing-package-s-parameter-estimation-for-diffeq-models/44082/16 "2020-08-05T15:37:32Z")

</div>

I think it’s a simple tilda, (shift+`)

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [August 5, 2020, 3:40pm UTC](https://discourse.julialang.org/t/turing-package-s-parameter-estimation-for-diffeq-models/44082/17 "2020-08-05T15:40:02Z")

</div>

I don’t understand quite what you mean. There is a LaTeX symbol `\tilde`, but if I type this + TAB, I get a different symbol, and it doesn’t work with `Turing`.

OK – I found it on the keyboard. But in general, it is quite “dangerous” to use very special keyboard symbols, because they may or may not exist on non-English keyboards.

---

<div class="post-metadata">

**Author:** ![arnaudmgh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arnaudmgh/32/7134_2.png) [@arnaudmgh](https://discourse.julialang.org/u/arnaudmgh)\
**Post date:** [August 5, 2020, 4:54pm UTC](https://discourse.julialang.org/t/turing-package-s-parameter-estimation-for-diffeq-models/44082/19 "2020-08-05T16:54:40Z")

</div>

if I remember well alt+N on several non-querty keyboards. Yes, keyboard layouts are a problem in programming, I am using US-QUERTY now, much easier than AZERTY for programming in my experience.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [August 5, 2020, 7:37pm UTC](https://discourse.julialang.org/t/turing-package-s-parameter-estimation-for-diffeq-models/44082/20 "2020-08-05T19:37:15Z")

</div>

Here is something that puzzles me… I ran both the “Fitting Lotka-Volterra with DiffEqBayes” (solid red + solid blue) procedure and the Turing model formulation (dotted green + dotted orange) twice. To save time, I only ran `1_000` samples for both methods. The resulting parameter distributions are very different:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/b/1/b1d622d6ee77f0112e2a5dc001fbe6a86d1b041c.png)

I have also done this with `10_000` samples – the results may differ substantially each time I run the Bayes fit.

_Questions_:

1. Do I do something wrong?
2. If I don’t do anything wrong, I don’t see how one can trust the results???

Here is my code… Method 1:

```julia
for i in 1:2
    bayes_est = turing_inference(prob,Tsit5(),tvec,data_lv,priors,num_samples=1_000)
    density!(fg_a,bayes_est.value.data[:,13,1],lc=LC[i])
    density!(fg_b,bayes_est.value.data[:,14,1],lc=LC[i])
    density!(fg_g,bayes_est.value.data[:,15,1],lc=LC[i])
    density!(fg_d,bayes_est.value.data[:,16,1],lc=LC[i])
    density!(fg_s,bayes_est.value.data[:,17,1],lc=LC[i])
end

```

…and Method 2:

```julia
for i in 1:2
    chain = sample(model, NUTS(.65), 1_000)
    density!(fg_a,chain.value.data[:,13,1],lc=LC[i],ls=:dot)
    density!(fg_b,chain.value.data[:,14,1],lc=LC[i],ls=:dot)
    density!(fg_g,chain.value.data[:,15,1],lc=LC[i],ls=:dot)
    density!(fg_d,chain.value.data[:,16,1],lc=LC[i],ls=:dot)
    density!(fg_s,chain.value.data[:,17,1],lc=LC[i],ls=:dot)
end

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 5, 2020, 9:35pm UTC](https://discourse.julialang.org/t/turing-package-s-parameter-estimation-for-diffeq-models/44082/21 "2020-08-05T21:35:35Z")

</div>

Looks like it had a bad initial chain one of the two times.

[Next page](https://discourse.julialang.org/t/turing-package-s-parameter-estimation-for-diffeq-models/44082.md?page=2)