Calculate Uncertainty for fitted parameter from least squares fit

How can I get the uncertainty for each of the fitted parameters after doing a least squares curve fit? I used tools-basic fitting- quadratic, but I could do the fit using lsqcurvefit or some other function if that is easier.

Each of the data points that were used for the curve fit had standard error associated with them, but I think that I can somehow calculate the uncertainty for each fitted parameter based off the residuals of the fit. Is this correct?

[Bernoulli Lizard]

I remember using Excel to calculate the chi squared values based off of the residuals, which somehow allowed to calculate the uncertainty for each parameter, one at a time. Is there any easier way to do this? Surely there must be a MATLAB function or routine for this by now...

Answer 1

[Matt J]

If your curve fit is unconstrained and your residual has uniform variance s2, then a common approximation to the covariance matrix of the parameters is

 Cov=inv(J'*J)*s2

where J is the Jacobian of the residual at the solution. Both LSQCURVEFIT and LSQNONLIN return the Jacobian as an optional output argument.

[Bernoulli Lizard]

Thank you, that's certainly a start.

But what does the covariance tell me about the uncertainty? Ultimately I want to calculate the uncertainty for a y value from the equation, given any x (so that I can draw error bars on discrete y values along the fitted line). I thought the way to do this would be to calculate the uncertainty on each fitted parameter and use propagation of uncertainty to calculate the uncertainty on y. Does this make sense?

thanks

[Bernoulli Lizard]

Also, I do not understand what s2 is here.

[Matt J]

No, it doesn't make sense. The curve parameter estimates are derived from y, so the uncertainty in y propagates into the parameters, not the other way around. If you don't know the uncertainty in y to begin with, you have no basis upon which to figure out the uncertainty in the parameters.

s2 is the variance of the errors in y(i). In my initial post above, I suggested that you might also view this as the variance of your residuals. This assumes, of course, that your curve fit is pretty close to the true y(i).

[Bernoulli Lizard]

So since s2 is the variance of the errors in y(i) (assuming that the curve fit is close, which I think it is), then sqrt(s2) is the uncertainty in y(i), right?

How can I actually calculate s2?

[Bernoulli Lizard]

Is the sqrt(norm of residuals) what I am looking for?

And should the uncertainty of the fit be constant, for every x value? The error for the original data points was obviously not constant, although I did not use a weighted least squares fit...

[Matt J]

Only you know what you are looking for. The term "uncertainty" is vague and without a precise technical definition.

What I can tell you is that if the errors in y(i) are Gaussian, with standard deviation sqrt(s2), then there is a 95% certainty that y(i) will be a distance of 2*sqrt(s2) from their mean. The mean of y(i) is presumably the fitted curve, or pretty close.

[Bernoulli Lizard]

I realize that 'uncertainty' is a very vague term. I think that what you told me to calculate is what I want to use, though.

But just to make sure I am interpreting this correctly: Is what MATLAB gives as 'resnorm', or the 'norm of residuals' what you are calling s2?

Thanks

[Matt J]

No. resnorm is something else.

[Bernoulli Lizard]

ok, so what is s2?

[Matt J]

I guess you're looking for me to rephrase what I told you (twice) before, but I'm not sure in what way: s2 is the variance of your residuals.

[Bernoulli Lizard]

I understand what it is mathematically, I just don't understand how to get that number from MATLAB. I've read the help documentation and I can't find out how to just get a list of all of the residuals so that I can find the variance of them.

[Matt J]

How about

 s2=var(y-curve);

where y(i) and curve(i) are the data points y and the fitted curve respectively?

[Graham Baker]

Can you explain why it is necessary to multiply by s2? In the definition of the covariance matrix that I'm familiar with, it would simply be calculated as cov=inv(J'*J).

[Matt J]

Well, the covariance of the parameter estimates has to depend on the statistical variability of the curve data y somehow. inv(J'*J) alone has no dependence on y whatsoever.