# Compute standard errors of nonlinear regression parameters with maximum likelihood method

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?

• The usual way is from the Hessian. What's the conditional distribution (given $t$) of the angle in your model? -- what you wrote is only the expectation of $\theta(t)$, which is only half a model. Is the variability about the mean constant (constant absolute error) or related to the mean (different at different "true" angles)? – Glen_b -Reinstate Monica Jun 13 '17 at 1:07
• Did you try bootstrapping? – Yannis Vassiliadis Jun 13 '17 at 7:17
• @Glen_b I suppose that all errors are distributed normally, so with a costant variability about the mean (if I understood correctly, I'm actually at the first year of a physics course). YannisVassiliadis How can I apply bootstrapping to this particular sample? These are not repeated measures – Francesco Manzali Jun 13 '17 at 10:06
• Then any decent nonlinear least squares regression software should estimate approximate standard errors for you. ... ctd – Glen_b -Reinstate Monica Jun 13 '17 at 10:44
• ctd... If you must do it yourself, note that you can turn your model into a partially linear one by writing the sin term as a sum of sin and cos terms; where conditional on the exponential term that's actually linear, and given that, the exponential one should just be an optimization over a single parameter I think. i.e you can look at a grid of values for b, estimate the remaining parameters via linear least squares for each one, get the residual sum of squares & then locate the optimum for b (perhaps using something like Brent's method to hone in on the optimum b pretty accurately) – Glen_b -Reinstate Monica Jun 13 '17 at 10:49

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

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!

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!

• Thanks, that's something I can easily implement! One more question: since I am resampling the same initial sample, is the final standard deviation distorted? I guess that if I did this procedure with real data in each sample (effectively repeating the same experiment $K$ times), with measures actually different between samples, I would get a standard deviation bigger than that I get from bootstrapping. If so, how can I correct for this? – Francesco Manzali Jun 14 '17 at 13:49
• Not necessarily. As $K\to\infty$ then the 2 distributions should converge. Also, if you had $K$ samples of $N$ observations each, you would probably run the estimation procedure on all $K\times N$ observations. – Yannis Vassiliadis Jun 14 '17 at 20:35