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

**URL:** <https://discourse.julialang.org/t/2d-lotka-volterra-parameter-estimation-with-only-1d-data/13838>\
**Category:** New to Julia\
**Tags:** diffeq, optim\
**Created:** [August 21, 2018, 4:23pm UTC](https://discourse.julialang.org/t/2d-lotka-volterra-parameter-estimation-with-only-1d-data/13838 "2018-08-21T16:23:13Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![gregory](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gregory/32/18969_2.png) [@gregory](https://discourse.julialang.org/u/gregory)\
**Post date:** [August 21, 2018, 4:23pm UTC](https://discourse.julialang.org/t/2d-lotka-volterra-parameter-estimation-with-only-1d-data/13838/1 "2018-08-21T16:23:13Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 23, 2018, 3:26pm UTC](https://discourse.julialang.org/t/2d-lotka-volterra-parameter-estimation-with-only-1d-data/13838/2 "2018-08-23T15:26:43Z")

</div>

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

```julia
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]

```

---

<div class="post-metadata">

**Author:** ![Matthew\_Overlin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/matthew_overlin/32/10533_2.png) [@Matthew\_Overlin](https://discourse.julialang.org/u/Matthew_Overlin)\
**Post date:** [January 10, 2020, 9:12pm UTC](https://discourse.julialang.org/t/2d-lotka-volterra-parameter-estimation-with-only-1d-data/13838/3 "2020-01-10T21:12:54Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 10, 2020, 10:07pm UTC](https://discourse.julialang.org/t/2d-lotka-volterra-parameter-estimation-with-only-1d-data/13838/4 "2020-01-10T22:07:42Z")

</div>

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

```julia
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))

```

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [January 21, 2020, 4:38pm UTC](https://discourse.julialang.org/t/2d-lotka-volterra-parameter-estimation-with-only-1d-data/13838/5 "2020-01-21T16:38:12Z")

</div>

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.

 ![Screenshot 2020-01-21 at 17.29.15](https://global.discourse-cdn.com/julialang/original/3X/1/2/12f69dd0677793289a323a2c2607aec7619084db.png)  
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.)

---

<div class="post-metadata">

**Author:** ![Matthew\_Overlin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/matthew_overlin/32/10533_2.png) [@Matthew\_Overlin](https://discourse.julialang.org/u/Matthew_Overlin)\
**Post date:** [January 22, 2020, 5:22am UTC](https://discourse.julialang.org/t/2d-lotka-volterra-parameter-estimation-with-only-1d-data/13838/6 "2020-01-22T05:22:59Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [February 1, 2020, 10:23am UTC](https://discourse.julialang.org/t/2d-lotka-volterra-parameter-estimation-with-only-1d-data/13838/7 "2020-02-01T10:23:04Z")

</div>

I put the code on the master branch of our package [BridgeSDEInference](https://github.com/mmider/BridgeSDEInference.jl)

The inference script is at  
[https://github.com/mmider/BridgeSDEInference.jl/blob/master/scripts/inference/lotka\_volterra.jl](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.

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [February 6, 2020, 2:28pm UTC](https://discourse.julialang.org/t/2d-lotka-volterra-parameter-estimation-with-only-1d-data/13838/8 "2020-02-06T14:28:37Z")

</div>

@Matthew_Overlin I forgot to tag you.