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

I want to fit some parameters and some initial conditions of an ODE system to data.
I am using Turing based on this tutorial.

Example code:

``````@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:

``````@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?

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