2D Lotka Volterra parameter estimation with only 1D data

gregory · August 21, 2018, 4:23pm

I am testing this with a toy problem first using a simple 2D Lotka Volterra. I feed in the actual params and then make some pseudo time series data. I then want to pretend I only have knowledge of 1 population and feed this as my objective function. Is there built-in functionality for this? If not, how could I go about putting this in.

tic()
using BlackBoxOptim
using DifferentialEquations
using RecursiveArrayTools # for VectorOfArray
using Optim
f1 = @ode_def_nohes LotkaVolterraTest begin
dx = x*(1 - x - Ay)
dy = rhoy*(1 - B*x - y)
end A B rho

u0 = [1.0;1.0]
tspan = (0.0,10.0)
p = [0.2,0.5,0.3]
prob = ODEProblem(f1,u0,tspan,p)
sol = solve(prob,Tsit5())
t = collect(linspace(0,10,200))
randomized = VectorOfArray([(sol(t[i]) + .01randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)

Here I would like to convert data → data[1,:] and then build a cost function to measure against this 1D time series, but plugging this into build_loss_objective returns an error.

cost_function = build_loss_objective(prob,Tsit5(),L2Loss(t,data),
maxiters=10000,verbose=false)

bound1 = Tuple{Float64, Float64}[(0,3),(0,3),(0,3)]
result = bboptimize(cost_function;SearchRange = bound1, MaxSteps = 5e4)

I would also like to do this with BFGS using a local optimizer:

tic()
using DifferentialEquations
using RecursiveArrayTools # for VectorOfArray
using Optim
f1 = @ode_def_nohes LotkaVolterraTest begin
dx = x*(1 - x - Ay)
dy = rhoy*(1 - B*x - y)
end A B rho

u0 = [1.0;1.0]
tspan = (0.0,10.0)
p = [0.2,0.5,0.3]
prob = ODEProblem(f1,u0,tspan,p)
sol = solve(prob,Tsit5())
t = collect(linspace(0,10,200))
randomized = VectorOfArray([(sol(t[i]) + .01randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)
cost_function = build_loss_objective(prob,Tsit5(),L2Loss(t,data),
maxiters=10000,verbose=false)

lower = [0.0,0.0,0.0]
upper = [3.0,3.0,3.0]
result_bfgs = optimize(cost_function, [1.5,1.7,1.6], lower, upper, Fminbox{BFGS}())

println(result_bfgs)
println(result_bfgs.minimizer)
toc()

ChrisRackauckas · August 23, 2018, 3:26pm

You can use save_idxs to make saveat only save the indices you want, and have that match up with the data. Example:

using DifferentialEquations, DiffEqParamEstim
using RecursiveArrayTools # for VectorOfArray
using Optim
f1 = @ode_def_nohes LotkaVolterraTest begin
  dx = x*(1 - x - A*y)
  dy = rho*y*(1 - B*x - y)
end A B rho

u0 = [1.0;1.0]
tspan = (0.0,10.0)
p = [0.2,0.5,0.3]
prob = ODEProblem(f1,u0,tspan,p)
sol = solve(prob,Tsit5())
t = collect(range(0,stop=10,length=200))
data = reshape([(sol(t[i];idxs=1) + .01randn()) for i in 1:length(t)],1,200)

cost_function = build_loss_objective(prob,Tsit5(),L2Loss(t,data),save_idxs = [1],
                                    maxiters=10000,verbose=false)
result_bfgs = optimize(cost_function, [.5,.7,.6], BFGS())

println(result_bfgs)
println(result_bfgs.minimizer) # [0.203781, 0.525022, 0.293092]


# Works with a patch on v1.0 to be released, slightly more efficient
cost_function = build_loss_objective(prob,Tsit5(),L2Loss(t,data),save_idxs = 1,
                                    maxiters=10000,verbose=false)
result_bfgs = optimize(cost_function, [.5,.7,.6], BFGS())

println(result_bfgs)
println(result_bfgs.minimizer) # [0.203781, 0.525022, 0.293092]

Matthew_Overlin · January 10, 2020, 9:12pm

Chris, when running your example, I get this error:

ERROR: LoadError: LoadError: UndefVarError: @ode_def_nohes not defined

I copied and pasted your example. Is there something else I should have done.
Thanks for the example.
Matt

ChrisRackauckas · January 10, 2020, 10:07pm

This example uses an old DSL that we don’t recommend any more. I would instead write the differential equation as:

function f(du,u,p,t)
  du[1] = p[1]*u[1] - p[2]*u[1]*u[2] #prey
  du[2] = -p[3]*u[2] + p[4]*u[1]*u[2] #predator
end
u0 = [1.0;1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1.0]
prob = ODEProblem(f,u0,tspan,p) 
t = collect(range(0, stop=10, length=200))

mschauer · January 21, 2020, 4:38pm

Hi @gregory, @Matthew_Overlin, is this still a hot topic for you?
We are having some success with our work on parameter inference of partially (1d data) observed multivariate ODE/SDE models and are thinking of using the Lotka-Volterra example for a tutorial/example.

The picture shows 50 very noisy 1D observations of a hidden Lotka-Volterra trajectory (black) with (α, β, γ, δ) = (1.5,1.0, 3.0,1.0) with additive noise (SDE). The red dots are the unknown true location of the process at observation times and the vertical red bars show the corresponding noisy 1D observations (vertical bars as the vertical position is unknown). The length of the dashed lines therefore indicates the observations errors.

This would be a setting where we get some reasonable estimates (α, β, γ, δ) = (1.59 0.84 (3.0) 0.95) (not all parameters are identifiable from partial observations so we fix γ to 3.)

Matthew_Overlin · January 22, 2020, 5:22am

I am still very interested in this topic. I would love to learn more about your example that you have here. Why isn’t one of the parameters possible to estimate? I could look through existing literature if you have good recommendations.

Also, I would love to try what you have, if you are ready and willing to share, to see how it might work for my purposes.

You have my attention.
-Matt

mschauer · February 1, 2020, 10:23am

I put the code on the master branch of our package BridgeSDEInference

The inference script is at
https://github.com/mmider/BridgeSDEInference.jl/blob/master/scripts/inference/lotka_volterra.jl

We are working on documentation and tutorials (and the paper!), but I commented the code and this should already give some idea.

mschauer · February 6, 2020, 2:28pm

@Matthew_Overlin I forgot to tag you.