Monte Carlo Parameter Estimation From Data

Chris Rackauckas

First you want to create a problem which solves multiple problems at the same time. This is the Monte Carlo Problem. When the parameter estimation tools say it will take any DEProblem, it really means ANY DEProblem!

So, let's get a Monte Carlo problem setup that solves with 10 different initial conditions.

using DifferentialEquations, DiffEqParamEstim, Plots, Optim

# Monte Carlo Problem Set Up for solving set of ODEs with different initial conditions

# Set up Lotka-Volterra system
function pf_func(du,u,p,t)
  du[1] = p[1] * u[1] - p[2] * u[1]*u[2]
  du[2] = -3 * u[2] + u[1]*u[2]
end
p = [1.5,1.0]
prob = ODEProblem(pf_func,[1.0,1.0],(0.0,10.0),p)
ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true
timespan: (0.0, 10.0)
u0: [1.0, 1.0]

Now for a MonteCarloProblem we have to take this problem and tell it what to do N times via the prob_func. So let's generate N=10 different initial conditions, and tell it to run the same problem but with these 10 different initial conditions each time:

# Setting up to solve the problem N times (for the N different initial conditions)
N = 10;
initial_conditions = [[1.0,1.0], [1.0,1.5], [1.5,1.0], [1.5,1.5], [0.5,1.0], [1.0,0.5], [0.5,0.5], [2.0,1.0], [1.0,2.0], [2.0,2.0]]
function prob_func(prob,i,repeat)
  ODEProblem(prob.f,initial_conditions[i],prob.tspan,prob.p)
end
monte_prob = MonteCarloProblem(prob,prob_func=prob_func)
MonteCarloProblem with problem ODEProblem

We can check this does what we want by solving it:

# Check above does what we want
sim = solve(monte_prob,Tsit5(),num_monte=N)
plot(sim)

nummonte=N means "run N times", and each time it runs the problem returned by the probfunc, which is always the same problem but with the ith initial condition.

Now let's generate a dataset from that. Let's get data points at every t=0.1 using saveat, and then convert the solution into an array.

# Generate a dataset from these runs
data_times = 0.0:0.1:10.0
sim = solve(monte_prob,Tsit5(),num_monte=N,saveat=data_times)
data = Array(sim)
2×101×10 Array{Float64,3}:
[:, :, 1] =
 1.0  1.06108   1.14403   1.24917   1.37764   …  0.956979  0.983561  1.0337
6 
 1.0  0.821084  0.679053  0.566893  0.478813     1.35559   1.10629   0.9063
71

[:, :, 2] =
 1.0  1.01413  1.05394  1.11711   …  1.05324  1.01309  1.00811  1.03162
 1.5  1.22868  1.00919  0.833191     2.08023  1.70818  1.39973  1.14803

[:, :, 3] =
 1.5  1.58801   1.70188   1.84193   2.00901   …  2.0153    2.21084   2.4358
9 
 1.0  0.864317  0.754624  0.667265  0.599149     0.600942  0.549793  0.5136
79

[:, :, 4] =
 1.5  1.51612  1.5621   1.63555   1.73531   …  1.83823   1.98545   2.15958 
 1.5  1.29176  1.11592  0.969809  0.850159     0.771089  0.691421  0.630025

[:, :, 5] =
 0.5  0.531705  0.576474  0.634384  0.706139  …  9.05366   9.4006   8.83911
 1.0  0.77995   0.610654  0.480565  0.380645     0.809382  1.51708  2.82619

[:, :, 6] =
 1.0  1.11027   1.24238   1.39866   1.58195   …  0.753108  0.748814  0.7682
84
 0.5  0.411557  0.342883  0.289812  0.249142     1.73879   1.38829   1.1093
2 

[:, :, 7] =
 0.5  0.555757  0.623692  0.705084  0.80158   …  8.11216   9.10671   9.9217
 
 0.5  0.390449  0.30679   0.24286   0.193966     0.261298  0.455937  0.8788
1

[:, :, 8] =
 2.0  2.11239   2.24921   2.41003   2.59433   …  3.223     3.47362   3.7301
4 
 1.0  0.909749  0.838025  0.783532  0.745339     0.739471  0.765597  0.8130
86

[:, :, 9] =
 1.0  0.969326  0.971358  1.00017  …  1.25065  1.1012   1.01733  0.979306
 2.0  1.63445   1.33389   1.09031     3.02671  2.52063  2.07502  1.69807 

[:, :, 10] =
 2.0  1.92148  1.88215  1.87711  1.90264  …  2.15079   2.27938   2.43105
 2.0  1.80195  1.61405  1.4426   1.2907      0.957221  0.884827  0.82948

Here, data[i,j,k] is the same as sim[i,j,k] which is the same as sim[k]i,j. So data[i,j,k] is the jth timepoint of the ith variable in the kth trajectory.

Now let's build a loss function. A loss function is some loss(sol) that spits out a scalar for how far from optimal we are. In the documentation I show that we normally do loss = L2Loss(t,data), but we can bootstrap off of this. Instead lets build an array of N loss functions, each one with the correct piece of data.

# Building a loss function
losses = [L2Loss(data_times,data[:,:,i]) for i in 1:N]
10-element Array{DiffEqParamEstim.L2Loss{StepRangeLen{Float64,Base.TwicePre
cision{Float64},Base.TwicePrecision{Float64}},Array{Float64,2},Nothing,Noth
ing,Nothing},1}:
 DiffEqParamEstim.L2Loss{StepRangeLen{Float64,Base.TwicePrecision{Float64},
Base.TwicePrecision{Float64}},Array{Float64,2},Nothing,Nothing,Nothing}(0.0
:0.1:10.0, [1.0 1.06108 … 0.983561 1.03376; 1.0 0.821084 … 1.10629 0.906371
], nothing, nothing, nothing, nothing)
 DiffEqParamEstim.L2Loss{StepRangeLen{Float64,Base.TwicePrecision{Float64},
Base.TwicePrecision{Float64}},Array{Float64,2},Nothing,Nothing,Nothing}(0.0
:0.1:10.0, [1.0 1.01413 … 1.00811 1.03162; 1.5 1.22868 … 1.39973 1.14803], 
nothing, nothing, nothing, nothing)   
 DiffEqParamEstim.L2Loss{StepRangeLen{Float64,Base.TwicePrecision{Float64},
Base.TwicePrecision{Float64}},Array{Float64,2},Nothing,Nothing,Nothing}(0.0
:0.1:10.0, [1.5 1.58801 … 2.21084 2.43589; 1.0 0.864317 … 0.549793 0.513679
], nothing, nothing, nothing, nothing)
 DiffEqParamEstim.L2Loss{StepRangeLen{Float64,Base.TwicePrecision{Float64},
Base.TwicePrecision{Float64}},Array{Float64,2},Nothing,Nothing,Nothing}(0.0
:0.1:10.0, [1.5 1.51612 … 1.98545 2.15958; 1.5 1.29176 … 0.691421 0.630025]
, nothing, nothing, nothing, nothing) 
 DiffEqParamEstim.L2Loss{StepRangeLen{Float64,Base.TwicePrecision{Float64},
Base.TwicePrecision{Float64}},Array{Float64,2},Nothing,Nothing,Nothing}(0.0
:0.1:10.0, [0.5 0.531705 … 9.4006 8.83911; 1.0 0.77995 … 1.51708 2.82619], 
nothing, nothing, nothing, nothing)   
 DiffEqParamEstim.L2Loss{StepRangeLen{Float64,Base.TwicePrecision{Float64},
Base.TwicePrecision{Float64}},Array{Float64,2},Nothing,Nothing,Nothing}(0.0
:0.1:10.0, [1.0 1.11027 … 0.748814 0.768284; 0.5 0.411557 … 1.38829 1.10932
], nothing, nothing, nothing, nothing)
 DiffEqParamEstim.L2Loss{StepRangeLen{Float64,Base.TwicePrecision{Float64},
Base.TwicePrecision{Float64}},Array{Float64,2},Nothing,Nothing,Nothing}(0.0
:0.1:10.0, [0.5 0.555757 … 9.10671 9.9217; 0.5 0.390449 … 0.455937 0.87881]
, nothing, nothing, nothing, nothing) 
 DiffEqParamEstim.L2Loss{StepRangeLen{Float64,Base.TwicePrecision{Float64},
Base.TwicePrecision{Float64}},Array{Float64,2},Nothing,Nothing,Nothing}(0.0
:0.1:10.0, [2.0 2.11239 … 3.47362 3.73014; 1.0 0.909749 … 0.765597 0.813086
], nothing, nothing, nothing, nothing)
 DiffEqParamEstim.L2Loss{StepRangeLen{Float64,Base.TwicePrecision{Float64},
Base.TwicePrecision{Float64}},Array{Float64,2},Nothing,Nothing,Nothing}(0.0
:0.1:10.0, [1.0 0.969326 … 1.01733 0.979306; 2.0 1.63445 … 2.07502 1.69807]
, nothing, nothing, nothing, nothing) 
 DiffEqParamEstim.L2Loss{StepRangeLen{Float64,Base.TwicePrecision{Float64},
Base.TwicePrecision{Float64}},Array{Float64,2},Nothing,Nothing,Nothing}(0.0
:0.1:10.0, [2.0 1.92148 … 2.27938 2.43105; 2.0 1.80195 … 0.884827 0.82948],
 nothing, nothing, nothing, nothing)

So losses[i] is a function which computes the loss of a solution against the data of the ith trajectory. So to build our true loss function, we sum the losses:

loss(sim) = sum(losses[i](sim[i]) for i in 1:N)
loss (generic function with 1 method)

As a double check, make sure that loss(sim) outputs zero (since we generated the data from sim). Now we generate data with other parameters:

prob = ODEProblem(pf_func,[1.0,1.0],(0.0,10.0),[1.2,0.8])
function prob_func(prob,i,repeat)
  ODEProblem(prob.f,initial_conditions[i],prob.tspan,prob.p)
end
monte_prob = MonteCarloProblem(prob,prob_func=prob_func)
sim = solve(monte_prob,Tsit5(),num_monte=N,saveat=data_times)
loss(sim)
10108.695792027418

and get a non-zero loss. So we now have our problem, our data, and our loss function... we have what we need.

Put this into buildlossobjective.

obj = build_loss_objective(monte_prob,Tsit5(),loss,num_monte=N,
                           saveat=data_times)
(::DiffEqParamEstim.DiffEqObjective{getfield(DiffEqParamEstim, Symbol("##29
#34")){Nothing,Bool,Int64,typeof(DiffEqParamEstim.STANDARD_PROB_GENERATOR),
Base.Iterators.Pairs{Symbol,Any,Tuple{Symbol,Symbol},NamedTuple{(:num_monte
, :saveat),Tuple{Int64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Ba
se.TwicePrecision{Float64}}}}},DiffEqBase.MonteCarloProblem{DiffEqBase.ODEP
roblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},DiffEq
Base.ODEFunction{true,typeof(Main.WeaveSandBox26.pf_func),LinearAlgebra.Uni
formScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,N
othing,Nothing},Nothing,DiffEqBase.StandardODEProblem},typeof(Main.WeaveSan
dBox26.prob_func),getfield(DiffEqBase, Symbol("##282#288")),getfield(DiffEq
Base, Symbol("##284#290")),Array{Any,1}},OrdinaryDiffEq.Tsit5,typeof(Main.W
eaveSandBox26.loss),Nothing},getfield(DiffEqParamEstim, Symbol("##33#39"))}
) (generic function with 2 methods)

Notice that I added the kwargs for solve into this. They get passed to an internal solve command, so then the loss is computed on N trajectories at data_times.

Thus we take this objective function over to any optimization package. I like to do quick things in Optim.jl. Here, since the Lotka-Volterra equation requires positive parameters, I use Fminbox to make sure the parameters stay positive. I start the optimization with [1.3,0.9], and Optim spits out that the true parameters are:

lower = zeros(2)
upper = fill(2.0,2)
result = optimize(obj, lower, upper, [1.3,0.9], Fminbox(BFGS()))
Results of Optimization Algorithm
 * Algorithm: Fminbox with BFGS
 * Starting Point: [1.3,0.9]
 * Minimizer: [1.5000000000581937,1.0000000001633538]
 * Minimum: 7.119929e-16
 * Iterations: 4
 * Convergence: true
   * |x - x'| ≤ 0.0e+00: true 
     |x - x'| = 0.00e+00 
   * |f(x) - f(x')| ≤ 0.0e+00 |f(x)|: true
     |f(x) - f(x')| = 0.00e+00 |f(x)|
   * |g(x)| ≤ 1.0e-08: false 
     |g(x)| = 1.07e-06 
   * Stopped by an increasing objective: true
   * Reached Maximum Number of Iterations: false
 * Objective Calls: 193
 * Gradient Calls: 193
result
Results of Optimization Algorithm
 * Algorithm: Fminbox with BFGS
 * Starting Point: [1.3,0.9]
 * Minimizer: [1.5000000000581937,1.0000000001633538]
 * Minimum: 7.119929e-16
 * Iterations: 4
 * Convergence: true
   * |x - x'| ≤ 0.0e+00: true 
     |x - x'| = 0.00e+00 
   * |f(x) - f(x')| ≤ 0.0e+00 |f(x)|: true
     |f(x) - f(x')| = 0.00e+00 |f(x)|
   * |g(x)| ≤ 1.0e-08: false 
     |g(x)| = 1.07e-06 
   * Stopped by an increasing objective: true
   * Reached Maximum Number of Iterations: false
 * Objective Calls: 193
 * Gradient Calls: 193

Optim finds one but not the other parameter.

I would run a test on synthetic data for your problem before using it on real data. Maybe play around with different optimization packages, or add regularization. You may also want to decrease the tolerance of the ODE solvers via

obj = build_loss_objective(monte_prob,Tsit5(),loss,num_monte=N,
                           abstol=1e-8,reltol=1e-8,
                           saveat=data_times)
result = optimize(obj, lower, upper, [1.3,0.9], Fminbox(BFGS()))
Results of Optimization Algorithm
 * Algorithm: Fminbox with BFGS
 * Starting Point: [1.3,0.9]
 * Minimizer: [1.500743476201907,1.001238477622136]
 * Minimum: 4.163900e-02
 * Iterations: 5
 * Convergence: true
   * |x - x'| ≤ 0.0e+00: true 
     |x - x'| = 0.00e+00 
   * |f(x) - f(x')| ≤ 0.0e+00 |f(x)|: true
     |f(x) - f(x')| = 0.00e+00 |f(x)|
   * |g(x)| ≤ 1.0e-08: false 
     |g(x)| = 1.07e-06 
   * Stopped by an increasing objective: false
   * Reached Maximum Number of Iterations: false
 * Objective Calls: 224
 * Gradient Calls: 224
result
Results of Optimization Algorithm
 * Algorithm: Fminbox with BFGS
 * Starting Point: [1.3,0.9]
 * Minimizer: [1.500743476201907,1.001238477622136]
 * Minimum: 4.163900e-02
 * Iterations: 5
 * Convergence: true
   * |x - x'| ≤ 0.0e+00: true 
     |x - x'| = 0.00e+00 
   * |f(x) - f(x')| ≤ 0.0e+00 |f(x)|: true
     |f(x) - f(x')| = 0.00e+00 |f(x)|
   * |g(x)| ≤ 1.0e-08: false 
     |g(x)| = 1.07e-06 
   * Stopped by an increasing objective: false
   * Reached Maximum Number of Iterations: false
 * Objective Calls: 224
 * Gradient Calls: 224

if you suspect error is the problem. However, if you're having problems it's most likely not the ODE solver tolerance and mostly because parameter inference is a very hard optimization problem.

Appendix

This tutorial is part of the DiffEqTutorials.jl repository, found at: https://github.com/JuliaDiffEq/DiffEqTutorials.jl

To locally run this tutorial, do the following commands:

using DiffEqTutorials
DiffEqTutorials.weave_file("ode_extras","04-monte_carlo_parameter_estim.jmd")

Computer Information:

Julia Version 1.1.1
Commit 55e36cc308 (2019-05-16 04:10 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-3770 CPU @ 3.40GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, ivybridge)

Package Information:

Status `~/.julia/environments/v1.1/Project.toml`
[7e558dbc-694d-5a72-987c-6f4ebed21442] ArbNumerics 0.5.4
[6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf] BenchmarkTools 0.4.2
[be33ccc6-a3ff-5ff2-a52e-74243cff1e17] CUDAnative 2.2.0
[3a865a2d-5b23-5a0f-bc46-62713ec82fae] CuArrays 1.0.2
[55939f99-70c6-5e9b-8bb0-5071ed7d61fd] DecFP 0.4.8
[abce61dc-4473-55a0-ba07-351d65e31d42] Decimals 0.4.0
[ebbdde9d-f333-5424-9be2-dbf1e9acfb5e] DiffEqBayes 1.1.0
[eb300fae-53e8-50a0-950c-e21f52c2b7e0] DiffEqBiological 3.8.2
[459566f4-90b8-5000-8ac3-15dfb0a30def] DiffEqCallbacks 2.5.2
[f3b72e0c-5b89-59e1-b016-84e28bfd966d] DiffEqDevTools 2.9.0
[1130ab10-4a5a-5621-a13d-e4788d82bd4c] DiffEqParamEstim 1.6.0
[055956cb-9e8b-5191-98cc-73ae4a59e68a] DiffEqPhysics 3.1.0
[6d1b261a-3be8-11e9-3f2f-0b112a9a8436] DiffEqTutorials 0.1.0
[0c46a032-eb83-5123-abaf-570d42b7fbaa] DifferentialEquations 6.4.0
[31c24e10-a181-5473-b8eb-7969acd0382f] Distributions 0.20.0
[497a8b3b-efae-58df-a0af-a86822472b78] DoubleFloats 0.9.1
[f6369f11-7733-5829-9624-2563aa707210] ForwardDiff 0.10.3
[c91e804a-d5a3-530f-b6f0-dfbca275c004] Gadfly 1.0.1
[7073ff75-c697-5162-941a-fcdaad2a7d2a] IJulia 1.18.1
[4138dd39-2aa7-5051-a626-17a0bb65d9c8] JLD 0.9.1
[23fbe1c1-3f47-55db-b15f-69d7ec21a316] Latexify 0.8.2
[eff96d63-e80a-5855-80a2-b1b0885c5ab7] Measurements 2.0.0
[961ee093-0014-501f-94e3-6117800e7a78] ModelingToolkit 0.2.0
[76087f3c-5699-56af-9a33-bf431cd00edd] NLopt 0.5.1
[2774e3e8-f4cf-5e23-947b-6d7e65073b56] NLsolve 4.0.0
[429524aa-4258-5aef-a3af-852621145aeb] Optim 0.18.1
[1dea7af3-3e70-54e6-95c3-0bf5283fa5ed] OrdinaryDiffEq 5.8.1
[65888b18-ceab-5e60-b2b9-181511a3b968] ParameterizedFunctions 4.1.1
[91a5bcdd-55d7-5caf-9e0b-520d859cae80] Plots 0.25.1
[d330b81b-6aea-500a-939a-2ce795aea3ee] PyPlot 2.8.1
[731186ca-8d62-57ce-b412-fbd966d074cd] RecursiveArrayTools 0.20.0
[90137ffa-7385-5640-81b9-e52037218182] StaticArrays 0.11.0
[f3b207a7-027a-5e70-b257-86293d7955fd] StatsPlots 0.11.0
[c3572dad-4567-51f8-b174-8c6c989267f4] Sundials 3.6.1
[1986cc42-f94f-5a68-af5c-568840ba703d] Unitful 0.15.0
[44d3d7a6-8a23-5bf8-98c5-b353f8df5ec9] Weave 0.9.0
[b77e0a4c-d291-57a0-90e8-8db25a27a240] InteractiveUtils
[37e2e46d-f89d-539d-b4ee-838fcccc9c8e] LinearAlgebra
[44cfe95a-1eb2-52ea-b672-e2afdf69b78f] Pkg