[SciPy-Dev] On the leastsq/curve_fit method

Gianfranco Durin g.durin@inrim...
Mon Oct 3 10:06:11 CDT 2011

<snip>
> The minimization of $\chi^{2}_{}$above is sometimes called weighted
> least squares in which case the inverse quantities 1/ e 2 are called
> the weights. Clearly this is simply a different word for the same
> thing, but in practice the use of these words sometimes means that
> the interpretation of e 2 as variances or squared errors is not
> straightforward. The word weight often implies that only the
> relative weights are known (point two is twice as important as
> point one'') in which case there is apparently an unknown overall
> normalization factor. Unfortunately the parameter errors coming out
> of such a fit will be proportional to this factor, and the user must
> be aware of this in the formulation of his problem.
> '''
> (I don't quite understand the last sentence.)
>
> M = σ^2 W^(-1), where σ^2 is estimated by residual sum of squares
> from weighted regression.
>
> W only specifies relative errors, the assumption is that the
> covariance matrix of the errors is *proportional* to W. The scaling
> is arbitrary. If the scale of W changes, then the estimated residual
> sum of squares from weighted regression will compensate for it. So,
> rescaling W has no effect on the covariance of the parameter
> estimates.
>
> I checked in Greene: Econometric Analysis, and briefly looked at the
> SAS description. It looks like weighted least squares is always with
> automatic scaling, W is defined only as relative weights.
>
> All I seem to be able to find is weighted least squares with
> automatic scaling (except for maybe some two-step estimators).
>
> Curve_fit calculates
>
> M = W \sigma^2 W^(-1) = \sigma^2
>
>
> If I remember correctly (calculated from the transformed model) it
> should be:
>
> the cov of the parameter estimates is s^2 (X'WX)^(-1)
> error estimates should be s^2 * W
>
> where W = diag(1/curvefit_sigma**2) unfortunate terminology for
> curve_fit's sigma or intentional ? (as I mentioned in the other
>
> Josef
>
> _______________________________________________
>
> Gianfranco,
> Can you please provide some Python and R (or SAS) code to show what
> you mean?
> ...
> Bruce

Bruce, Josef and the others,
after reading a few books, searching over many website, and tried many different software packages (and, with the help of many colleagues), I came to a conclusion which should fix the question:

Weights and data errors (or uncertainties) CAN be different concepts, as written above. It depends on the user...
In particular:
1. If he/she thinks the sigma JUST as weights and NOT as standard-deviation of ydata, cov_x MUST be multiplied by the residual standard deviation
2. If the user thinks the sigma as standard-deviation of ydata (i.e. measurement errors), which are by the way ALSO good to weight the data themself, then cov_x DOES NOT NEED to be multiplied by the residual standard deviation.

I tried a few packages, and found they assume one of the two options by default without 'asking' (or making aware of) the user. In particular (check using the very same data and errors):
Option 1: SAS, R, gnuplot, octave...
Option 2: Profit, Origin, ...

And mathematica? In the HOWTO: Fit Models with Measurement Errors (see below), mathematica makes the difference between weights and measurement errors, so the user can decide how to use his/her sigma.

I think we should make this distinction explicit also in our curve_fit, and report it on the leastsq doc.

Gianfranco

============================================================
from Mathematica:
"Particularly in the physical sciences,it is common to use measurement errors as weights to incorporate measured variation into the fitting. Weights have a relative effect on the parameter estimates, but an error variance still needs to be estimated in weighted regression, and this impacts error estimates for results."

1. When using Weights alone, the variance scale is estimated using the default method [i.e. our s_sq]. Error estimates will depend on both the weights and the estimated variance scale. However, if the weights are from measurement errors, you would want error estimates to depend solely on the weights.
It is important to note that weights do not change the fitting or error estimates. For example, multiplying all weights by a constant increases the estimated variance,but does not change the parameter estimates or standard errors. (Ps. This is what I meant saying that the parameters' errors
....

2. For measurements errors, you want standard errors to be computed only from the weights.... While weights have an impact on parameter estimates, the variance estimate itself does not."