0
$\begingroup$

I'm examining a code in C++ for a nonlinear fit. It is basically a Levenberg Marquardt routine you can find on Netlib or elsewhere. The last step is estimating the errors of the parameters that are fitted. From literature, I know that the variance of my parameters, $a_j$, is given by $\sigma^2 a_j = C_{j,j}$ where $C$ is the inverse of the Hessian. In the code this is exactly what happens, i.e. by a Gauss-Jordan algorithm the Hessian is inverted. What puzzles me is the last step to calculate the relative error:

$\delta a_j = \sqrt{\frac{f(x, \vec{a}_\text{fit})}{50 - N} \vert C_{j,j} \vert } \frac{1}{a_j}$

The value $f$ is the minimum value, i.e. the minimum sum of squares. $N$ is the number of parameters that are to be fitted. $\vert C_{j,j} \vert$ is the absolute entry in the inverted Hessian. $a_j$ is the parameter for which the error is to be estimated.

From a dimensional point of view, the result is a percentage, which is fine, but the denominator $50 - N$ puzzles me. At the end, there's a test whether $\delta a_j$ is larger than 50 %, but I can't see that this is connected with the shown calculation. If someone has an idea of the denominator, thanks for a hint.

Improve this question
$\endgroup$
3
  • $\begingroup$ Can you cite &/or link to where you got this formula? Is it what is coded into the C++ function you're using? $\endgroup$
    – gung - Reinstate Monica
    Apr 18, 2015 at 17:22
  • $\begingroup$ As @Brian says in his answer, 50 must equal the sample size, and the sample size minus the number of fit parameters (N, here), is the degrees of freedom which appears in the denominator. $\endgroup$
    – Harvey Motulsky
    Apr 18, 2015 at 22:29
  • $\begingroup$ Yes, @gung, it is a C++ function I was examining; the explanation of Brian is correct. I wasn't realizing this; even though the number of data points is accessible at that point, 50 was chosen as a fixed number. $\endgroup$
    – Clemens
    Apr 19, 2015 at 6:49

1 Answer 1

Reset to default
2
$\begingroup$

Assuming that $f(x,a_{fit})$ is the sum of squares of the residuals, then

$s=\sqrt{\frac{f(x,a_{fit})}{50-N}}$

is an estimate of the standard deviation of the measured data points. The denominator of $50-N$ takes into account that $f$ is the sum of the squares of 50 residuals. You're assuming here that each of the 50 measurements has normally distributed errors that are independent of each other wish standard deviation $s$. The $-N$ in the denominator accounts approximately for the fact that you've fit $N$ parameters to the data. If you didn't include the $-N$, then you'd be effectively over fitting the data and underestimating $s$.

This is discussed in many books that cover nonlinear regression. A fairly elementary treatment is in the book by Bevington and Robinson.

You should be aware that there are many approximations involved with this approach to estimating a confidence interval for the fitted parameter. It's easy to find nonlinear regression problems for which those approximations don't work well. You'd be well advised to to test these approximations before you use them. Bevington and Robinson discuss many of the problems that can occur.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Let's assume one wants to fit a Gaussian to a grayscale image - i.e. having lots of data points and only very few parameters - would you say that this is still a "sensible" approach, or is one better off by just taking $\sqrt{|C_{jj}|}$ (notation from above) for an error estimate? $\endgroup$
    – Clemens
    Apr 25, 2015 at 19:05
  • $\begingroup$ No, I wouldn't just use $\sqrt{C_{j,j}}$. $\endgroup$
    – Brian Borchers
    Apr 25, 2015 at 19:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.