 # 2D Lotka Volterra parameter estimation with only 1D data

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 = 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(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 = 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(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()

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 = ,
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]
``````

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

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

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 = p*u - p*u*u #prey
du = -p*u + p*u*u #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))
``````