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

Specific Domains Modelling & Simulations
, ,
1

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?

2

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