Hello all,

We need to fit a quadratic to a set of data, and use the coefficients and their standard deviations. Using the "lfit" routine, we can obtain the correct coefficients straightforwardly. However, we are a bit confused by the covariance matrix.

Our understanding is that we get the standard deviations by taking the square root of the values on the diagonal of the covariance matrix. But, for example, the following simple data set fitted to 2 coefficients:

x = 1, y = 3
x = 2, y = 5
x = 3, y = 7
x = 4, y = 9
x = 5, y = 11

returns the coefficients 1 and 2 as expected, but the covariance matrix is

1.1 -0.3
-0.3 0.1

the diagonal of which we would expect to be 0 for both coefficients, as is the case for the straight-line routine "fit".

We would be very grateful for your assistance.

Many thanks

Edward Wiles
Tessella Support Services, UK
edward.wiles@tessella.com

Hello again,

It appears that the trick is to run "lfit" twice.

1) Run "lfit" on the data, with all values of the parameter "sig" set to a default 1.

2) Use the returned coefficients to calculate the residual standard deviation of the data in a straightforward manner. (Ignore the covariance matrix.) For a straight-line fit, this residual standard deviation is exactly the parameter "chi2" in the routine "fit".

3) Run "lfit" again, this time with all values of "sig" set to the residual standard deviation calculated in 2) above.

4) Find the square root of each value on the diagonal of the covariance matrix returned in 3) above. These are the standard deviations of the coefficients. [Note that the coefficients returned in 3) are the same as those returned in 1).]

We'd be grateful for any further thoughts or corrections to the above.

Many thanks

Edward