[SciPy-User] Unexpected covariance matrix from scipy.optimize.curve_fit
josef.pktd at gmail.com
josef.pktd at gmail.com
Wed Aug 31 12:10:52 EDT 2011
On Wed, Aug 31, 2011 at 11:09 AM, Christoph Deil
<deil.christoph at googlemail.com> wrote:
>
> On Aug 30, 2011, at 11:25 PM, josef.pktd at gmail.com wrote:
>
> On Tue, Aug 30, 2011 at 4:15 PM, Christoph Deil
> <deil.christoph at googlemail.com> wrote:
>
> I noticed that scipy.optimize.curve_fit returns parameter errors that don't
>
> scale with sigma, the standard deviation of ydata, as I expected.
>
> Here is a code snippet to illustrate my point, which fits a straight line to
>
> five data points:
>
> import numpy as np
>
> from scipy.optimize import curve_fit
>
> x = np.arange(5)
>
> y = np.array([1, -2, 1, -2, 1])
>
> sigma = np.array([1, 2, 1, 2, 1])
>
> def f(x, a, b):
>
> return a + b * x
>
> popt, pcov = curve_fit(f, x, y, p0=(0.42, 0.42), sigma=sigma)
>
> perr = np.sqrt(pcov.diagonal())
>
> print('*** sigma = {0} ***'.format(sigma))
>
> print('popt: {0}'.format(popt))
>
> print('perr: {0}'.format(perr))
>
> I get the following result:
>
> *** sigma = [1 2 1 2 1] ***
>
> popt: [ 5.71428536e-01 1.19956213e-08]
>
> perr: [ 0.93867933 0.40391117]
>
> Increasing sigma by a factor of 10,
>
> sigma = 10 * np.array([1, 2, 1, 2, 1])
>
> I get the following result:
>
> *** sigma = [10 20 10 20 10] ***
>
> popt: [ 5.71428580e-01 -2.27625699e-09]
>
> perr: [ 0.93895295 0.37079075]
>
> The best-fit values stayed the same as expected.
>
> But the error on the slope b decreased by 8% (the error on the offset a
>
> didn't change much)
>
> I would have expected fit parameter errors to increase with increasing
>
> errors on the data!?
>
> Is this a bug?
>
> No bug in the formulas. I tested all of them when curve_fit was added.
>
> However in your example the numerical cov lacks quite a bit of
> precision. Trying your example with different starting values, I get a
> 0.05 difference in your perr (std of parameter estimates).
>
> Trying smaller xtol and ftol doesn't change anything. (?)
>
> Making ftol = 1e-15 very small I get a different wrong result:
> popt: [ 5.71428580e-01 -2.27625699e-09]
> perr: [ 0.92582011 0.59868281]
> What do I have to do to get a correct answer (say to 5 significant digits)
> from curve_fit for this simple example?
>
> Since it's linear
>
> import scikits.statsmodels.api as sm
>
> x = np.arange(5.)
>
> y = np.array([1, -2, 1, -2, 1.])
>
> sigma = np.array([1, 2, 1, 2, 1.])
>
> res = sm.WLS(y, sm.add_constant(x, prepend=True), weights=1./sigma**2).fit()
>
> res.params
>
> array([ 5.71428571e-01, 1.11022302e-16])
>
> res.bse
>
> array([ 0.98609784, 0.38892223])
>
> res = sm.WLS(y, sm.add_constant(x, prepend=True),
> weights=1./(sigma*10)**2).fit()
>
> res.params
>
> array([ 5.71428571e-01, 1.94289029e-16])
>
> res.bse
>
> array([ 0.98609784, 0.38892223])
>
> rescaling doesn't change parameter estimates nor perr
>
> This is what I don't understand.
> Why don't the parameter estimate errors increase with increasing errors
> sigma on the data points?
> If I have less precise measurements, the model parameters should be less
> constrained?!
> I was using MINUIT before I learned Scipy and the error definition for a
> chi2 fit given in the MINUIT User Guide
> http://wwwasdoc.web.cern.ch/wwwasdoc/minuit/node7.html
> as well as the example results here
> http://code.google.com/p/pyminuit/wiki/GettingStartedGuide
> don't mention the factor s_sq that is used in curve_fit to scale pcov.
> Is the error definition in the MINUIT manual wrong?
> Can you point me to a web resource that explains why the s_sq factor needs
> to be applied to the covariance matrix?
It's standard text book information, but Wikipedia seems to be lacking
a bit in this.
for the linear case
http://en.wikipedia.org/wiki/Ordinary_least_squares#Assuming_normality
cov_params = sigma^2 (X'X)^{-1}
for the non-linear case with leastsq, X is replaced by Jacobian,
otherwise everything is the same.
However, in your minuit links I saw only the Hessian mentioned (from
very fast skimming the pages)
With maximum likelihood, the inverse Hessian is the complete
covariance matrix, no additional multiplication is necessary.
Essentially, these are implementation details depending on how the
estimation is calculated, and there are various ways of numerically
approximating the Hessian.
That's why this is described for optimize.leastsq (incorrectly as
Chuck pointed out) and but not in optimize.curve_fit.
With leastsquares are maximum likelihood, rescaling both y and
f(x,params) has no effect on the parameter estimates, it's just like
changing units of y, meters instead of centimeters.
I guess scipy.odr would work differently, since it is splitting up the
errors between y and x's, but I never looked at the details.
>
> Josef
>
>
>
> PS: I've attached a script to fit the two examples using statsmodels, scipy
> and minuit (applying the s_sq factor myself).
> Here are the results I get (who's right for the first example? why does
> statsmodels only return on parameter value and error?):
> """Example from
> http://code.google.com/p/pyminuit/wiki/GettingStartedGuide"""
> x = np.array([1 , 2 , 3 , 4 ])
> y = np.array([1.1, 2.1, 2.4, 4.3])
> sigma = np.array([0.1, 0.1, 0.2, 0.1])
> statsmodels.api.WLS
> popt: [ 1.04516129]
> perr: [ 0.0467711]
> scipy.optimize.curve_fit
> popt: [ 8.53964011e-08 1.04516128e+00]
> perr: [ 0.27452122 0.09784324]
that's what I get with example 1 when I run your script,
I don't know why you have one params in your case
(full_output threw an exception in curve_fit with scipy.__version__ '0.9.0'
statsmodels.api.WLS
popt: [ -6.66133815e-16 1.04516129e+00]
perr: [ 0.33828314 0.12647671]
scipy.optimize.curve_fit
popt: [ 8.53964011e-08 1.04516128e+00]
perr: [ 0.27452122 0.09784324]
> minuit
> popt: [-4.851674617611934e-14, 1.0451612903225629]
> perr: [ 0.33828315 0.12647671]
statsmodels and minuit agree pretty well
> """Example from
> http://mail.scipy.org/pipermail/scipy-user/2011-August/030412.html"""
> x = np.arange(5)
> y = np.array([1, -2, 1, -2, 1])
> sigma = 10 * np.array([1, 2, 1, 2, 1])
> statsmodels.api.WLS
> popt: [ 5.71428571e-01 7.63278329e-17]
> perr: [ 0.98609784 0.38892223]
> scipy.optimize.curve_fit
> popt: [ 5.71428662e-01 -8.73679511e-08]
> perr: [ 0.97804034 0.3818681 ]
> minuit
> popt: [0.5714285714294132, 2.1449508835758024e-13]
> perr: [ 0.98609784 0.38892223]
statsmodels and minuit agree,
my guess is that the jacobian calculation of leastsq (curve_fit) is
not very good in these examples. Maybe trying Dfun or the other
options, epsfcn, will help.
I was trying to see whether I get better results calculation the
numerical derivatives in a different way, but had to spend the time
fixing bugs.
(NonlinearLS didn't work correctly with weights.)
Josef
>
>
>
>
> Looking at the source code I see that scipy.optimize.curve_fit multiplies
>
> the pcov obtained from scipy.optimize.leastsq by a factor s_sq:
>
> https://github.com/scipy/scipy/blob/master/scipy/optimize/minpack.py#L438
>
> if (len(ydata) > len(p0)) and pcov is not None:
>
> s_sq = (func(popt, *args)**2).sum()/(len(ydata)-len(p0))
>
> pcov = pcov * s_sq
>
> If so is it possible to add an explanation to
>
> http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html
>
> that pcov is multiplied with this s_sq factor and why that will give correct
>
> errors?
>
> After I noticed this issue I saw that this s_sq factor is mentioned in the
>
> cov_x return parameter description of leastsq,
>
> but I think it should be explained in curve_fit where it is applied, maybe
>
> leaving a reference in the cov_x leastsq description.
>
> Also it would be nice to mention the full_output option in the curve_fit
>
> docu, I only realized after looking at the source code that this was
>
> possible.
>
> Christoph
>
> _______________________________________________
>
> SciPy-User mailing list
>
> SciPy-User at scipy.org
>
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>
More information about the SciPy-User
mailing list