# Fitting parameters together with initial conditions of an ODE to real data

**URL:** <https://discourse.julialang.org/t/fitting-parameters-together-with-initial-conditions-of-an-ode-to-real-data/93255>\
**Category:** Modelling & Simulations\
**Tags:** question, turing, differentialequation\
**Created:** [January 20, 2023, 9:52am UTC](https://discourse.julialang.org/t/fitting-parameters-together-with-initial-conditions-of-an-ode-to-real-data/93255 "2023-01-20T09:52:26Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![moeinkh88](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moeinkh88/32/29607_2.png) [@moeinkh88](https://discourse.julialang.org/u/moeinkh88)\
**Post date:** [January 20, 2023, 9:52am UTC](https://discourse.julialang.org/t/fitting-parameters-together-with-initial-conditions-of-an-ode-to-real-data/93255/1 "2023-01-20T09:52:26Z")

</div>

I want to fit some parameters and some initial conditions of an ODE system to data.  
I am using Turing based on this [tutorial](https://turing.ml/dev/tutorials/10-bayesian-differential-equations/).

Example code:

```julia
@model function fitprob(data)
    # Prior distributions.
    σ ~ InverseGamma(2, 3)
	# μ ~ truncated(Normal(0, 1); lower=0, upper=1/665205)
    Λ ~ truncated(Normal(μ* 665205, 665205); lower=μ* 665205, upper=665205)
	μp ~ truncated(Normal(0, 1); lower=0, upper=1)
	# ϕ1 ~ truncated(Normal(0, 1); lower=0, upper=1)
	ϕ2 ~ truncated(Normal(0, 1); lower=0.002, upper=.75)
	β1 ~ truncated(Normal(0, 1); lower=0, upper=1)
	β2 ~ truncated(Normal(0, 1); lower=0, upper=1)
	δ ~ truncated(Normal(0, 1); lower=0, upper=1)
	ψ ~ truncated(Normal(0, 1); lower=0, upper=1)
	# ω ~ truncated(Normal(0, 1); lower=0, upper=1)
	σ2 ~ truncated(Normal(0, 1); lower=0.006, upper=0.01)
	# γS ~ truncated(Normal(0, 1); lower=0, upper=1)
	# γA ~ truncated(Normal(0, 1); lower=0, upper=1)
	ηS ~ truncated(Normal(0, .5); lower=0, upper=.5)
	ηA ~ truncated(Normal(0, .5); lower=0, upper=.5)
	IA ~ truncated(Normal(0,Λ/μ - 665205); lower=0, upper=min(Λ/μ - 665205,65188))
	P ~ truncated(Normal(0,Λ*(ηA+ηS)/(μ*μp)); lower=0, upper=min(Λ*(ηA+ηS)/(μ*μp),65188))

    # Simulate model.
	p=[Λ,μ,μp,ϕ1,ϕ2,β1,β2,δ,ψ,ω,σ2,γS,γA,ηS,ηA]
	X0=[665188,0,IA,17,0,P,0,665188+17+IA]
	prob1 = ODEProblem(F, X0, tSpan, par)
    predict = solve(prob1,alg_hints=[:stiff]; p=p, saveat=1, save_idxs=[5,7])

    # Observations.
    for i in 1:length(predict)
        data[i, :] ~ MvNormal(predict[i], σ^2 * I)
    end

    return nothing
end

model = fitprob([TR TD]) # joint vector data 

# Sample 3 independent chains with forward-mode automatic differentiation (the default).
chain = sample(model, NUTS(0.65), MCMCSerial(), 5000, 3; progress=false)

```

I am wondering if it is OK to define the problem inside the model function since its initial condition is varying.

Is it possible to fix the problem and call it inside the `solve` by setting new initial conditions? For example something like:

```julia
@model function fitprob(data, prob)
[...]
predict = solve(prob,alg_hints=[:stiff]; p=p, u0=X0, saveat=1, save_idxs=[5,7])
[...]

```

Is there any better way of setting?

How about controlling the bound of parameters based on other parameters? Is this `truncated(Normal(0,Λ*(ηA+ηS)/(μ*μp)); lower=0, upper=min(Λ*(ηA+ηS)/(μ*μp),65188))` too dumb?

---

<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:** [January 20, 2023, 11:10am UTC](https://discourse.julialang.org/t/fitting-parameters-together-with-initial-conditions-of-an-ode-to-real-data/93255/2 "2023-01-20T11:10:06Z")

</div>

> [@moeinkh88](#):
>
> I am wondering if it is OK to define the problem inside the model function since its initial condition is varying.

Yes

> Is there any better way of setting?

You can also use `remake` and set parameters and initial conditions accordingly like `prob = remake(prob; p = p, u0 = X0)`

> How about controlling the bound of parameters based on other parameters? Is this `truncated(Normal(0,Λ*(ηA+ηS)/(μ*μp)); lower=0, upper=min(Λ*(ηA+ηS)/(μ*μp),65188))` too dumb?

I think that’s the way to do it