Compute standard errors of nonlinear regression parameters with maximum likelihood method

Question

I have a set of experimental data ${(t_i,\theta_i)}_{i=1\,\dots\,N}$ of the angular positions $\theta$ at the istants $t$ measured from a lightly damped harmonic oscillator, and i want fit them with the theoretical model: $$ \theta(t) = a\,e^{-bt}\,sin(ct+d) $$ I'm using the maximum likelihood method for parameter estimation, and I can compute the most likely parameters given a certain set of data. However, I can't find a closed solution for the problem of maximizing $\mathcal{L}$, so i cannot use the usual error propagation formula in order to find the standard errors of the fitted parameters.

What should I do to compute the parameters' errors?

Answer 1

I wrote a little Python helper to help with this problem (see here). You can use the fit.get_vcov() function to get the standard errors of the parameters. It uses automatic differentiation to compute the Hessian and uses that to compute the standard errors of the best-fit parameters.

Background

You can look through the slides here, but I will explain it as best as I can. You want the standard errors of the best-fit parameters, which is the same as the standard deviation of the best-fit parameters. 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))

$\sigma$ and $H$ can give you $\Sigma$. Just read the diagonals off this matrix for the standard errors!

Answer 2

So bootstrapping would be a good idea. The basic advantage is that you don't have to worry about closed form/analytical solutions, you can just resample. The basic algorithm would be as follows:

Assuming you have a dataset with N observations and a vector $\hat{\theta}$ of estimated parameters.

Step 1) Create a new sample $N_k, k = 1, \dots, K$ by sampling N observations with replacement from your original dataset. It NEEDS to be with replacement, otherwise you always end up with the sample and the same $\hat{\theta}$. As you said, each observation can appear multiple times, or not at all.

Step 2) Run the estimation procedure, and store vector $\hat{\theta_k}$.

Step 3) Repeat K times

Repeat a large number of times, I think $K = 1000$ would suffice (or more if you want). We do that because when the Law of Large Numbers kicks in, we have $\hat{\theta} \xrightarrow{d} \theta$. So, this will give us the sampling distribution of the test statistic.

The standard deviation of the sampling distribution of $\hat{\theta}$ is the standard error of $\hat{\theta}$. Similarly, if you want to get confidence intervals, determine what level of significance ($\alpha$) you want, and based on that remove the tails of the distribution.

All standard statistical software will do the entire process for you. If you want better theoretical understanding, the wikipedia entry on bootstrapping is actually pretty good.

Sorry for the long answer, but I hope it helps!