adnls

Using automatic differentiation as implemented in the autograd package for analyzing non-linear curve fitting. Presently finds the Jacobian, Hessian, variance covariance matrix of fitted parameters for user-defined non-linear model.

Quick start/Example

Example taken from Niclas Börlin's lecture slides, which also has a clear explanation of the concepts involved. Download adnls.py from this repo into your working directory. Then the model and residual function are defined.

from adnls import fit
import autograd.numpy as np

# The model used for nonlinear regression
def model(t, p):
    y = p[0] * np.exp(p[1] * t)
    return y

# The residual function used by the `fit` class
def res(p, extra_pars):
    return model(extra_pars["t"], p) - extra_pars["data"]

Then define the data to be stored in the extra_pars variable, which is a dictionary. Best fitting values are assumed to have been found by some other method.

extra_pars = {}
extra_pars["t"] = np.array([1, 2, 4])
extra_pars["data"] = np.array([3, 4, 6])
x_best = np.array((2.4893671, 0.2210597), dtype=np.float64)

Define the instance of the fit class and then call the requisite functions.

this_fit = fit(
    res_func=res,
    bestfit_pars=x_best,
    extra_pars=extra_pars,
)  # create instance of fit class
print("Hessian")
print(this_fit.get_H(x_best))
print(this_fit.get_H())
print("Jacobian")
print(this_fit.get_J(x_best))
print("Residual")
print(this_fit.residual())
print(this_fit.residual(x_best))
print("Vcov matrix")
print(this_fit.get_vcov())
print("Sd of best-fit paramters")
print(this_fit.get_sd_bf())

Output

Hessian
[[  9.83906461  74.29758072] 
 [ 74.29758072 651.85399352]]
[[  9.83906461  74.29758072] 
 [ 74.29758072 651.85399352]]
Jacobian
[[ 1.2473979   3.10523129]
 [ 1.55600152  7.74691796]
 [ 2.42114072 24.10843219]]
Residual
[ 0.10523129 -0.12654102  0.02710805]
[ 0.10523129 -0.12654102  0.02710805]
Vcov matrix
[[ 0.02029685 -0.00231341]
 [-0.00231341  0.00030636]]
Sd of best-fit paramters
[0.142467   0.01750314]

Background

We are interested finding the standard deviation (sd) of the best-fit parameters (or the standard errors of the best-fit parameters) and conducting significance tests on each parameter. The sd of the best fit parameters are given by the diagonal elements of the covariance matrix $\Sigma$. $\Sigma$ for non-linear regression is given by:

$$ \Sigma = \sigma^2 (H^{-1}) $$

where $\sigma$ is the standard deviation of the residuals and $H$ is the Hessian of the objective function (such as least squares or weighted least squares).

Finding standard deviation of the residuals, $\sigma$

If you don't know $\sigma$ from previous experiments, then you can estimate it as $\hat{\sigma}$ and use that estimated value to get $\Sigma = \hat{\sigma}^2 (H^{-1})$. It can be estimated with:

$$ \hat{\sigma} = \sqrt{\frac{f(x_{best})}{m-n}} $$

where $f(x_{best})$ is the best likelihood found by maximum-likelihood (aka best fit objective function). This can be something like the sum of squared residuals (SSE). $m$ is the number of parameters in your model. $n$ is the number of data points used to fit your model.

Finding the Hessian of the objective function, $H$

The Hessian is the same as the Jacobian of the gradient. I use the autograd.elementwise_grad and autograd.jacobian to compute $H$. See the snippet below:

from autograd import elementwise_grad as egrad
from autograd import jacobian
import autograd.numpy as np


def func(x):
    return np.sin(x[0]) * np.sin(x[1])


x_value = np.array([0.0, 0.0])  # has to be float, not int
H_f = jacobian(egrad(func))  # returns a function
print(H_f(x_value))

Together with $\sigma$ (or $\hat{\sigma}$), we can estimate $\Sigma$. From there, the standard errors of the paramters are just the diagonal elements of $\Sigma$.

$t$ value and testing for significance

The $t$ statistic for the $i^th$ is estimated as $x_{best,i}$ / $i^th$ element of diag($\Sigma$). This should be $t$ distributed with degrees of freedom = $m-n$.

Features

1. Compute Jacobian of the residuals with respect to fitted parameters using automatic differentiation.
2. Compute Hessian of the objective function with respect to fitted parameters using automatic differentiation.
3. Estimate standard deviation if it is not known.
4. Estimate variance covariance matrix of the fitted parmeters, derived the Hessian found by automatic differentiation.
5. Estimate the standard deviations of the fitted parameters, derived the Hessian found by automatic differentiation.

Requirements

Install the autograd package with pip install autograd. See here for further instructions.

Caveats

1. Model must use autograd.numpy and autograd.scipy where applicable.

2. Only the subset of numpy and scipy implemented by the autograd package will work. See here for list of things that will work.

3. Differential equation models don't appear to work at the moment. For example, consider the differential equation variant of the model above:

   from adnls import fit
   import autograd.numpy as np
   from autograd.scipy.integrate import odeint
   # from scipy.integrate import odeint
   
   
   def derivative(y, t, p):
       dy = p[0] * p[1] * np.exp(p[1] * t)
       return dy
   
   
   def model(t, p):
       y0 = np.array([0])
       sol = odeint(derivative, y0, t, args=(p,))
       print(sol)
       return sol[1:, 0] + p[0]
   
   
   def res(p, extra_pars):
       return model(extra_pars["t"], p) - extra_pars["data"]
       # [ 0.10523129 -0.12654102  0.02710805]
   
   
   if __name__ == "__main__":
       extra_pars = {}
       extra_pars["t"] = np.array([0, 1, 2, 4])
       extra_pars["data"] = np.array([3, 4, 6])
       x_best = np.array((2.4893671, 0.2210597), dtype=np.float64)
       print(res(x_best, extra_pars))
       this_fit = fit(
           res_func=res,
           bestfit_pars=x_best,
           extra_pars=extra_pars,
       )
       print("Jacobian")
       print(this_fit.get_J(x_best))

   This fails with

    ...
     File "C:\Users\username\Documents\projects\adnls\adnls.py", line 109, in <lambda>
       J = jacobian(lambda pars: self._res_func(pars, self._extra_pars))(x)
     File "diffeq.py", line 21, in res
       return model(extra_pars["t"], p) - extra_pars["data"]
     File "diffeq.py", line 15, in model
       sol = odeint(derivative, y0, t, args=(p,))
     File "C:\Users\username\Miniconda3\envs\ad\lib\site-packages\autograd\tracer.py", line 48, in f_wrapped
       return f_raw(*args, **kwargs)
     File "C:\Users\username\Miniconda3\envs\ad\lib\site-packages\scipy\integrate\odepack.py", line 244, in odeint
       int(bool(tfirst)))
   ValueError: setting an array element with a sequence.
   

   The model itself is fine, though. This can be verified by uncommeting # from scipy.integrate import odeint and commenting out from autograd.scipy.integrate import odeint. Not sure if it's a bug in my implementation.

Work-around

If you must use a differential equation model, considering using the up-to-date JAX library. Presently there is an experimental ODE solver in JAX at https://github.com/google/jax/blob/master/jax/experimental/ode.py .

• But this only works on Linux/macOS. See here for discussion on porting JAX for Windows - apparently this can be done with Windows Subsystem for Linux (WSL). (The challenge was getting jaxlib to compile on Windows)

To do

1. Try out WSL for differential equation models.
2. Add reference for background.