ODE parameter estimation on a subset of the parameters

General Usage
,
1

Hello, I’m a differential equation newbie.

I have an ODE model with 50 parameters and 10 compartments. I have data measured in only 3 compartments and I wish to fit 6 parameters. The 3 compartment problem is discussed here, but how do I fit my 6 parameters while holding the other 40 fixed?

For example, I modified the Lotka-Volterra model (docs) so the -3 becomes the second parameter of p:

function f(du,u,p,t)
  du[1] = dx = p[1]*u[1] - u[1]*u[2]
  du[2] = dy = p[2]*u[2] + u[1]*u[2]
end

#problem & initial condition
u0 = [1.0, 1.0]
tspan = (0.0, 10.0)
p = [1.5, -3]
prob = ODEProblem(f,u0,tspan,p)

#simulate data
t = collect(range(0,stop=10,length=200))
using RecursiveArrayTools # for VectorOfArray
randomized = VectorOfArray([(sol(t[i]) + .01randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)

#define cost function to calculate error in both compartments
cost_function = build_loss_objective(prob,Tsit5(),L2Loss(t, data),
                                     maxiters=10000,verbose=false, save_idxs=[1, 2])

function closure()
	# what goes here?
end

Question: How would I just fit the 1st parameter p[1] and hold p[2] constant at -3?

Possible method 1:
Following some recommendation on Slack, I am using save_idxs in build_loss_objective, and I’m trying to build a closure to only optimize the parameters I want. However, I can’t figure out what this closure is supposed to look like. Should I redefine a new ODE based on the subset of parameters I care about in this closure? (probably not).

Possible method 2:
The other way that will work is to redefine a new ODE model where only the parameters I want to estimate are parameters, and those that should be fixed are hard-coded into the equations. However, I have more than 50 parameters and I may wish to estimate different subsets of them in the future, so this method is probably not a good way either.

2

well, I have the following solution for a closure:

function closure(p)
	p[2] = -1.5
	err = cost_function(p)
	return err
end

result = optimize(cost_function, [1.0, -3.0])

julia> result.minimizer
2-element Array{Float64,1}:
  1.5000022369113704
 -2.9999430279237727

I suppose this works?

3

The easiest solution is just leave p2 as a number and not a parameter.
Otherwise I have an example at blackboxoptim where I deal exactly with this.

https://github.com/robertfeldt/BlackBoxOptim.jl/blob/master/examples/ODEwithFixParams.jl

I’m out of the country without computer so I’m not much help, but this example I think address your issue-

I used Julia 1.04 to run the code.

Take care-

A