[SciPy-dev] least squares error
Ondrej Certik
ondrej at certik.cz
Mon Mar 9 16:53:48 EDT 2009
On Mon, Mar 9, 2009 at 7:27 AM, <josef.pktd at gmail.com> wrote:
> On Mon, Mar 9, 2009 at 8:16 AM, Ondrej Certik <ondrej at certik.cz> wrote:
>> Hi,
>>
>> I think each time I used leastsq, I also needed to calculate the
>> errors in the fitted parameters. I use this method, that takes the
>> output of leastsq and returns the parameters+errors.
>>
>> def calc_error(args):
>> p, cov, info, mesg, success = args
>> chisq=sum(info["fvec"]*info["fvec"])
>> dof=len(info["fvec"])-len(p)
>> sigma = array([sqrt(cov[i,i])*sqrt(chisq/dof) for i in range(len(p))])
>> return p, sigma
>>
>> let's integrate this with leastsq? E.g. add a new key in the info dict?
>>
>> Ondrej
>> _______________________________________________
>
> That's close to what the new scipy.optimize.curve_fit does, except it
> returns the correctly scaled covariance matrix of the parameter
> estimate
>
> popt, pcov = curve_fit(func, x, yn)
> psigma = np.sqrt(np.diag(pcov)) ) #standard deviation of parameter estimates
>
> Note: using chisq=sum(info["fvec"]*info["fvec"]) saves one function
> call compared to the curve_fit implementation.
>
> I would prefer that optimize.leastsq stays a low level wrapper, and
> the interpretation and additional statistics are produced in higher
> level functions, such as curve_fit.
Yes, but given the complexity of my code above (e.g. it's trivial), I
think it could also be added to leastsq to the info dict, because it
already returns some stat. data.
>
> The higher level functions can take care of calculating the correct
> covariance of the parameter estimates for different cases, e.g. when
> using weights as in curve_fit, or generalized least squares, or
> sandwich estimates for the covariance.
>
> What I would like to see additional in optimize.leastsq is the
> Jacobian directly, since this can be used to test additional
> restrictions on the parameter estimates, or derive additional
> statistics. I still haven't figured out whether it is possible to
> recover it.
I thought the Jacobian is in the info dict as well.
On Mon, Mar 9, 2009 at 1:47 PM, <josef.pktd at gmail.com> wrote:
> On Mon, Mar 9, 2009 at 3:37 PM, Nils Wagner
> <nwagner at iam.uni-stuttgart.de> wrote:
>> Hi,
>>
>> there is another issue concerning least squares
>>
>> http://projects.scipy.org/numpy/ticket/937
>>
>> Nils
>
> ticket 937 is for numpy.linalg.lstsq
> while I assume Ondrej meant scipy.optimize.leastsq
Yep.
Ondrej
More information about the SciPy-Dev
mailing list