[SciPy-user] incorrect variance in optimize.curvefit and leastsq

Robert Kern robert.kern@gmail....
Fri Feb 13 01:06:13 CST 2009


On Fri, Feb 13, 2009 at 00:08,  <josef.pktd@gmail.com> wrote:
> On Fri, Feb 13, 2009 at 12:24 AM, Travis E. Oliphant
> <oliphant@enthought.com> wrote:
>> josef.pktd@gmail.com wrote:
>>> On Thu, Feb 12, 2009 at 11:09 PM, Travis E. Oliphant
>>> <oliphant@enthought.com> wrote:
>>>
>>>> josef.pktd@gmail.com wrote:
>>>>
>>>>> I just saw the new optimize.curvefit which provides a wrapper around
>>>>> optimize.leastsq
>>>>>
>>>>>
>>>>>
>>> the standard deviation of the error can be calculated and the
>>> corrected (this is written for the use from outside of curvefit):
>>>
>>> yhat = func(x,popt[0], popt[1])   # get predicted observations
>>> SSE = np.sum((y-yhat)**2)
>>> sig2 = SSE/(len(y)-len(popt))
>>> ecov = sig2*pcov       # this is the variance-covariance  matrix of
>>> the parameter estimates
>>>
>>>
>>> inside curvefit, this work (before the return):
>>>
>>> err  = func(popt, *args)
>>> SSE = np.sum((err)**2)
>>> sig2 = SSE / (len(ydata) - len(popt))
>>> pcov = sig2 * pcov
>
> I think for the weighted least squares problem the weights should go
> into the SSE calculation. I didn't find directly the reference, but I
> am somewhat confident that this is correct, from the analogy to the
> transformed
> model ynew = y*weight where weight_i = 1/sigma_i in the linear case.
> But it's too late today to try to check this.
>
> SSE = np.sum((err*weight)**2)

Yes. Basically, you want a Chi-squared statistic for the residuals.

scipy.odr does this scaling, for example.

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
  -- Umberto Eco


More information about the SciPy-user mailing list