How to estimate confidence interval of a least-squares fit parameters by means of numerical Jacobian

I am working on a complicated data fitting algorithm in Matlab. I have a problem with properly estimating the confidence intervals of my fit. I will describe my procedure in some detail, give some of my thoughts on the problem and subsequently formulate my question more precisely.

Disclaimer: I am a physicist, not a statistician. Please respond in technical, but not discipline-specific language if possible.

What do I do?

I have a complicated data set and an even more complicated model for it which I don't want to describe here since that is not relevant. I do a least squares fit of the data set with respect to $P$ - a vector of fit parameters (I am not using maximum likelihood estimation yet). This part works very well and I've tested it extensively on known (both real and artificial) data-sets.

I subsequently aim to follow this procedure in order to estimate the prediction bounds on my fit parameters. From now on I will use the same variable nomenclature so I recommend reading the procedure now.

Problem

My $f(x,P)$ - the fit function is purely numerical, and so I use a numerical estimate of $J_{f}$ (using this code). $J_f$ is then a $\texttt{length}(P) \times N$ matrix - Number of data points $\times$ number of fit parameters evaluated at the best fit parameters - $P'$. $$J_{f} = \left[\begin{array}{c c c c}\frac{\partial f(P',x_1)}{\partial a_1} & \frac{\partial f(P',x_2)}{\partial a_1} & \dots & \frac{\partial f(P',x_N)}{\partial a_1}\\ \frac{\partial f(P',x_1)}{\partial a_2} & \frac{\partial f(P',x_2)}{\partial a_2} & \dots & \frac{\partial f(P',x_N)}{\partial a_2}\\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial f(P',x_1)}{\partial a_k} & \frac{\partial f(P',x_2)}{\partial a_k} & \dots & \frac{\partial f(P',x_N)}{\partial a_k}\end{array}\right]$$ I then estimate the Hessian as $$H \approx J_{f}^TJ_{f}$$ such that a matrix element of $H$ is $$H_{ij} = \sum_{n = 1}^{N}\left( \frac{\partial f(P',x_n)}{\partial a_i} \times \frac{\partial f(P',x_n)}{\partial a_j} \right)$$ (sidenote: Hessian obtained like that seems to agree with numerical Hessian obtained using this code).

$H$ is the observed Fisher information matrix (according to this). Moreover an inverse of Fisher information matrix is an estimator of the asymptotic covariance matrix so $$C \approx \sigma_r H^{-1}$$ where $\sigma_r$ is an unbiased variance of the residuals.

The standard errors of P are therefore $$P_{se} = \sqrt{\sigma_r\, \texttt{diag}(H^{-1})} \times \texttt{tinv}(1-0.05/2,v)$$ where the last term is $\approx$ 1.96 - Student's t inverse cumulative distribution function factor for 95% confidence interval, and $v$ is the number of degrees of freedom. This assumes a normal distribution of the errors in the fitted data.

Questions

• Is this procedure at all valid?
• When estimating the Hessian matrix should it in fact be $$H \approx \frac{1}{N}J_{f}^{T}J_{f}$$ such that $$H_{ij} = \,{\bf\frac{1}{N}}\;\sum_{n = 1}^{N}\left( \frac{\partial f(P',x_n)}{\partial a_i} \times \frac{\partial f(P',x_n)}{\partial a_j} \right)$$ ? This would mean that H is invariant with respect to $N$
• Instead of getting the Jacobian of $f(x,P)$ (which is a vector valued function) should I instead use Jacobian of the root-mean-square error ($\varepsilon$) of my fit?

Root-mean-square error function: $$\varepsilon(P,x,y)= \frac{1}{N}\left(\sum\limits_{i=1}^N (f(x,P)-y)^2\right)^{\frac{1}{2}}$$ Jacobian: $$J_{\varepsilon} = \left[\begin{array}{c}\frac{\partial \varepsilon(P',x,y)}{\partial a_1} \\ \frac{\partial \varepsilon(P',x,y)}{\partial a_2}\\ \vdots\end{array}\right]$$ Standard errors: $$H \approx (J_e^TJ_e)$$ etc.

If I am understanding the question properly, you are asking about the non-linear least squares (NLS) model. Sticking to your notation, the NLS model is: \begin{align} y_n &= f(x_n,P) + \epsilon_n \end{align} Various regularity conditions are required on $f$, $X$, and $\epsilon$. The parameters are $P$, a $k$-vector. They are estimated by solving: \begin{align} min_P \sum_{n=1}^N \left(y_n-f(x_n,P) \right)^2 \end{align} Assuming that you have found the global minimum and that it is interior to whatever the feasible set is for $P$ ($\mathbb{R}^k$, I guess?), the necessary FOC are: \begin{align} J_f\left(Y-f(X,P)\right)=0 \end{align} The estimator, $P'$, defined as the interior solution to the minimization problem above and therefore solving the first-order condition immediately above, is consistent and asymptotically normal. It has an asymptotic variance of $\sigma^2_rH^{-1}$. The variance of the error in the original model is assumed to be $Var(\epsilon_n)=\sigma^2_r$, and it may be consistently estimated as $\hat{\sigma}_r^2=\frac{1}{N-k}\sum_{n=1}^N\left(y_n-f(x_n,P')\right)^2$.

These are standard results for the non-linear least squares model. My favorite reference for them is Amemiya, T (1985) Advanced Econometrics, section 4.3.

The upshot of all this is that, if you want a 95% confidence interval for, say, the $3^{rd}$ element of $P$, you can use: \begin{align} P_{3}' \pm 1.96\sqrt{\hat{\sigma}_r^2\left(H^{-1}\right)_{33}} \end{align}

Notice that one of the regularity conditions is that all the variances of the $\epsilon$ are equal and that all of the covariances among elements of the $\epsilon$ are zero. The model in the link you provided contemplates that different $\epsilon_n$ have different variances (i.e. heteroskedasticity).

If that characterizes your application, then you need to modify your variance matrix. You can use a so-called sandwich estimator, like this: \begin{align} C &= \left( J_f^T J_f \right)^{-1}J_f^T \hat{\Sigma} J_f \left( J_f^T J_f \right)^{-1}\\ \hat{\Sigma} &= diag \left(\left(y_n-f(x_n,P') \right)^2 \right) \end{align} Then, the 95% confidence interval for the $3^{rd}$ element of $P$, would be: \begin{align} P_{3}' \pm 1.96\sqrt{C_{33}} \end{align} Further modifications would be required if there were to be correlations among the various elements of $\epsilon$.

• Thank you, that is very helpful. It does not actually answer all my questions though. Given all your assumptions about my problem (which are basically correct), can I use the $H≈J^T_fJ_f$ estimate for $H$? Would it make sense to use $H≈\frac1NJ^T_fJ_f$? – MarcinKonowalczyk Feb 28 '17 at 10:36
• Also; would this procedure still work if you did not assume that the global minimum has been reached? – MarcinKonowalczyk Feb 28 '17 at 10:42
• Yes, you should use $H=J^T_f J_f$. It would not make sense to divide by $N$. This would cause the variance matrix to converge to something positive-definite (usually). That makes no sense. The variance is going to go to zero as $N$ goes to infinity, assuming some regularity conditions. – Bill Feb 28 '17 at 18:58
• Basically, the answer to your other question is no. It does not matter if you have not obtained the global minimum. But the reason this is true might bother you. To get the estimator $P'$ to be consistent for the true value of the parameter, you need your objective function (times $1/N$) to converge to a function with a unique minimum on the interior of the feasible set of parameter values. So, as $N$ goes to infinity, there is only one solution to the first order condition, eventually. So, the whole local/global problem goes away in the limit. – Bill Feb 28 '17 at 19:34
• For myself, when I estimate a model with multiple local maxima, I spend some effort to try to find the global maximum (though this is difficult to ensure, of course). My intuition is that the global maximum will have better small sample properties. I don't know of a proof for this, though. – Bill Feb 28 '17 at 19:41