# Parameter Estimation

Parameter estimation for ODE models, also known as dynamic data analysis, is provided by the DiffEq suite. Note these require that the problem is defined using a ParameterizedFunction.

## Recommended Methods

The recommended method is to use `build_loss_objective`

with the optimizer of your choice. This method can thus be paired with global optimizers from packages like NLopt.jl which can be much less prone to finding local minima than local optimization methods. Also, it allows the user to define the cost function in the way they choose as a function `loss(sol)`

, and thus can fit using any cost function on the solution, making it applicable to fitting non-temporal data and other types of problems. Also, `build_loss_objective`

works for all of the `DEProblem`

types, allowing it to optimize parameters on ODEs, SDEs, DDEs, DAEs, etc.

However, this method requires repeated solution of the differential equation. If the data is temporal data, the most efficient method is the `two_stage_method`

which does not require repeated solutions but is not as accurate. Usage of the `two_stage_method`

should have a post-processing step which refines using a method like `build_loss_objective`

.

### two_stage_method

The two-stage method is a collocation method for estimating parameters without requiring repeated solving of the differential equation. It does so by determining a smoothed estimated trajectory of the data and optimizing the derivative function and the data's timepoints to match the derivatives of the smoothed trajectory. This method has less accuracy than other methods but is much faster, and is a good method to try first to get in the general "good parameter" region, to then finish using one of the other methods.

```
function two_stage_method(prob::DEProblem,tpoints,data;kernel= :Epanechnikov,
loss_func = L2DistLoss,mpg_autodiff = false,
verbose = false,verbose_steps = 100)
```

### build_loss_objective

`build_loss_objective`

builds an objective function to be used with Optim.jl and MathProgBase-associated solvers like NLopt.

```
function build_loss_objective(prob::DEProblem,alg,loss_func
regularization=nothing;
mpg_autodiff = false,
verbose_opt = false,
verbose_steps = 100,
prob_generator = problem_new_parameters,
kwargs...)
```

The first argument is the `DEProblem`

to solve, and next is the `alg`

to use. The `alg`

must match the problem type, which can be any `DEProblem`

(ODEs, SDEs, DAEs, DDEs, etc.). `regularization`

defaults to nothing which has no regulariztion function. One can also choose `verbose_opt`

and `verbose_steps`

, which, in the optimization routines, will print the steps and the values at the steps every `verbose_steps`

steps. `mpg_autodiff`

uses autodifferentiation to define the derivative for the MathProgBase solver. The extra keyword arguments are passed to the differential equation solver.

#### The Loss Function

`loss_func(sol)`

is a function which reduces the problem's solution to a scalar which the optimizer will try to minimize. While this is very flexible, two convenience routines are included for fitting to data:

```
L2Loss(t,data;weight=nothing)
CostVData(t,data;loss_func = L2Loss,weight=nothing)
```

where `t`

is the set of timepoints which the data is found at, and `data`

are the values that are known. `L2Loss`

is an optimized version of the L2-distance. In `CostVData`

, one can choose any loss function from LossFunctions.jl or use the default of an L2 loss. The `weight`

is a vector of weights for the loss function which must match the size of the data.

Note that minimization of a weighted `L2Loss`

is equivalent to maximum likelihood estimation of a heteroskedastic Normally distributed likelihood.

#### Note About Loss Functions

For parameter estimation problems, it's not uncommon for the optimizers to hit unstable regions of parameter space. This causes warnings that the solver exited early, and the built-in loss functions like `L2Loss`

and `CostVData`

automatically handle this. However, if using a user-supplied loss function, you should make sure it's robust to these issues. One common pattern is to apply infinite loss when the integration is not successful. Using the retcodes, this can be done via:

```
function my_loss_function(sol)
tot_loss = 0.0
if any((s.retcode != :Success for s in sol))
tot_loss = Inf
else
# calculation for the loss here
end
tot_loss
end
```

#### The Regularization Function

The regularization can be any function of `p`

, the parameter vector:

`regularization(p)`

The `Regularization`

helper function builds a regularization using a penalty function `penalty`

from PenaltyFunctions.jl:

`Regularization(λ,penalty=L2Penalty())`

The regularization defaults to L2 if no penalty function is specified. `λ`

is the weight parameter for the addition of the regularization term.

#### The Problem Generator Function

The argument `prob_generator`

allows one to specify a function for generating new problems from a given parameter set. By default, this just builds a new version of `f`

that inserts all of the parameters. For example, for ODEs this is given by the dispatch on `DiffEqBase.problem_new_parameters`

which does the following:

```
function problem_new_parameters(prob::ODEProblem,p)
f = (t,u,du) -> prob.f(t,u,p,du)
uEltype = eltype(p)
u0 = [uEltype(prob.u0[i]) for i in 1:length(prob.u0)]
tspan = (uEltype(prob.tspan[1]),uEltype(prob.tspan[2]))
ODEProblem(f,u0,tspan)
end
```

`f = (t,u,du) -> prob.f(t,u,p,du)`

creates a new version of `f`

that encloses the new parameters. The element types for `u0`

and `tspan`

are set to match the parameters. This is required to make autodifferentiation work. Then the new problem with these new values is returned.

One can use this to change the meaning of the parameters using this function. For example, if one instead wanted to optimize the initial conditions for a function without parameters, you could change this to:

```
function my_problem_new_parameters(prob::ODEProblem,p)
uEltype = eltype(p)
tspan = (uEltype(prob.tspan[1]),uEltype(prob.tspan[2]))
ODEProblem(prob.f,p,tspan)
end
```

which simply matches the type for time to `p`

(once again, for autodifferentiation) and uses `p`

as the initial condition in the initial value problem.

### build_lsoptim_objective

`build_lsoptim_objective`

builds an objective function to be used with LeastSquaresOptim.jl.

`build_lsoptim_objective(prob,tspan,t,data;prob_generator = problem_new_parameters,kwargs...)`

The arguments are the same as `build_loss_objective`

.

### lm_fit

`lm_fit`

is a function for fitting the parameters of an ODE using the Levenberg-Marquardt algorithm. This algorithm is really bad and thus not recommended since, for example, the Optim.jl algorithms on an L2 loss are more performant and robust. However, this is provided for completeness as most other differential equation libraries use an LM-based algorithm, so this allows one to test the increased effectiveness of not using LM.

`lm_fit(prob::DEProblem,tspan,t,data,p0;prob_generator = problem_new_parameters,kwargs...)`

The arguments are similar to before, but with `p0`

being the initial conditions for the parameters and the `kwargs`

as the args passed to the LsqFit `curve_fit`

function (which is used for the LM solver). This returns the fitted parameters.

## Local Optimization Examples

We choose to optimize the parameters on the Lotka-Volterra equation. We do so by defining the function as a ParmaeterizedFunction:

```
f = @ode_def_nohes LotkaVolterraTest begin
dx = a*x - b*x*y
dy = -c*y + d*x*y
end a=>1.5 b=1.0 c=3.0 d=1.0
u0 = [1.0;1.0]
tspan = (0.0,10.0)
prob = ODEProblem(f,u0,tspan)
```

Notice that since we only used `=>`

for `a`

, it's the only free parameter. We create data using the numerical result with `a=1.5`

:

```
sol = solve(prob,Tsit5())
t = collect(linspace(0,10,200))
using RecursiveArrayTools # for VectorOfArray
randomized = VectorOfArray([(sol(t[i]) + .01randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)
```

Here we used `VectorOfArray`

from RecursiveArrayTools.jl to turn the result of an ODE into a matrix.

If we plot the solution with the parameter at `a=1.42`

, we get the following:

Notice that after one period this solution begins to drift very far off: this problem is sensitive to the choice of `a`

.

To build the objective function for Optim.jl, we simply call the `build_loss_objective`

function:

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

Note that we set `maxiters`

in a way that causes the differential equation solvers to error more quickly when in bad regions of the parameter space, speeding up the process. Now this cost function can be used with Optim.jl in order to get the parameters. For example, we can use Brent's algorithm to search for the best solution on the interval `[0,10]`

by:

```
using Optim
result = optimize(cost_function, 0.0, 10.0)
```

This returns `result.minimizer[1]==1.5`

as the best parameter to match the data. When we plot the fitted equation on the data, we receive the following:

Thus we see that after fitting, the lines match up with the generated data and receive the right parameter value.

We can also use the multivariate optimization functions. For example, we can use the `BFGS`

algorithm to optimize the parameter starting at `a=1.42`

using:

`result = optimize(cost_function, [1.42], BFGS())`

Note that some of the algorithms may be sensitive to the initial condition. For more details on using Optim.jl, see the documentation for Optim.jl.

Lastly, we can use the same tools to estimate multiple parameters simultaneously. Let's use the Lotka-Volterra equation with all parameters free:

```
f2 = @ode_def_nohes LotkaVolterraAll begin
dx = a*x - b*x*y
dy = -c*y + d*x*y
end a=>1.5 b=>1.0 c=>3.0 d=>1.0
u0 = [1.0;1.0]
tspan = (0.0,10.0)
prob = ODEProblem(f2,u0,tspan)
```

To solve it using LeastSquaresOptim.jl, we use the `build_lsoptim_objective`

function:

`cost_function = build_lsoptim_objective(prob,Tsit5(),L2Loss(t,data))`

The result is a cost function which can be used with LeastSquaresOptim. For more details, consult the documentation for LeastSquaresOptim.jl:

```
x = [1.3,0.8,2.8,1.2]
res = optimize!(LeastSquaresProblem(x = x, f! = cost_function,
output_length = length(t)*length(prob.u0)),
LeastSquaresOptim.Dogleg(),LeastSquaresOptim.LSMR(),
ftol=1e-14,xtol=1e-15,iterations=100,grtol=1e-14)
```

We can see the results are:

```
println(res.minimizer)
Results of Optimization Algorithm
* Algorithm: Dogleg
* Minimizer: [1.4995074428834114,0.9996531871795851,3.001556360700904,1.0006272074128821]
* Sum of squares at Minimum: 0.035730
* Iterations: 63
* Convergence: true
* |x - x'| < 1.0e-15: true
* |f(x) - f(x')| / |f(x)| < 1.0e-14: false
* |g(x)| < 1.0e-14: false
* Function Calls: 64
* Gradient Calls: 9
* Multiplication Calls: 135
```

and thus this algorithm was able to correctly identify all four parameters.

## More Algorithms (Global Optimization) via MathProgBase Solvers

The `build_loss_objective`

function builds an objective function which is able to be used with MathProgBase-associated solvers. This includes packages like IPOPT, NLopt, MOSEK, etc. Building off of the previous example, we can build a cost function for the single parameter optimization problem like:

```
f = @ode_def_nohes LotkaVolterraTest begin
dx = a*x - b*x*y
dy = -c*y + d*x*y
end a=>1.5 b=1.0 c=3.0 d=1.0
u0 = [1.0;1.0]
tspan = (0.0,10.0)
prob = ODEProblem(f,u0,tspan)
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)
obj = build_loss_objective(prob,Tsit5(),L2Loss(t,data),maxiters=10000)
```

We can now use this `obj`

as the objective function with MathProgBase solvers. For our example, we will use NLopt. To use the local derivative-free Constrained Optimization BY Linear Approximations algorithm, we can simply do:

```
using NLopt
opt = Opt(:LN_COBYLA, 1)
min_objective!(opt, obj)
(minf,minx,ret) = NLopt.optimize(opt,[1.3])
```

This finds a minimum at `[1.49997]`

. For a modified evolutionary algorithm, we can use:

```
opt = Opt(:GN_ESCH, 1)
min_objective!(opt, obj.cost_function2)
lower_bounds!(opt,[0.0])
upper_bounds!(opt,[5.0])
xtol_rel!(opt,1e-3)
maxeval!(opt, 100000)
(minf,minx,ret) = NLopt.optimize(opt,[1.3])
```

We can even use things like the Improved Stochastic Ranking Evolution Strategy (and add constraints if needed). This is done via:

```
opt = Opt(:GN_ISRES, 1)
min_objective!(opt, obj.cost_function2)
lower_bounds!(opt,[-1.0])
upper_bounds!(opt,[5.0])
xtol_rel!(opt,1e-3)
maxeval!(opt, 100000)
(minf,minx,ret) = NLopt.optimize(opt,[0.2])
```

which is very robust to the initial condition. The fastest result comes from the following:

```
using NLopt
opt = Opt(:LN_BOBYQA, 1)
min_objective!(opt, obj)
(minf,minx,ret) = NLopt.optimize(opt,[1.3])
```

For more information, see the NLopt documentation for more details. And give IPOPT or MOSEK a try!