[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