[SciPy-User] optimize.leastsq error estimation
Matt Newville
newville@cars.uchicago....
Thu Mar 21 21:38:25 CDT 2013
Hi Joe,
On Thu, Mar 21, 2013 at 3:10 PM, Joe Philip Ninan <indiajoe@gmail.com> wrote:
> Hi,
> I want to confirm whether i understood the documentation of optimize.leastsq
> correctly.
> (
> http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.html
> )
> For estimating the error in the fitted parameter.
> The documentation says :
> "This matrix (cov_x) must be multiplied by the residual variance to get the
> covariance of the parameter estimates"
> i.e. Multiply the output cov_x with the "reduced chi square". Right? and
> take sqrt of corresponding diagonal element from covariance matrix.
> Where "residual chi square" is the sum of {[(fitted function -
> data)/sigma]**2} /N
> and where N is degrees of freedom [= len(data)-number of parameters]
> And sigma is the error in each data point.
>
> I mainly want to confirm whether the divide by sigma is required or not.
> (I couldn't find the definition of "residual variance" in documentation, so
> i assumed it is reduced chi square).
>
> Thanking you,
> Joe
I think you're correct, but there may be some confusion in
terminology. The standard errors are the square-root of the diagonal
elements of (covar * reduced_chi_square), where reduced_chi_square is
(((data-model)/uncertainty)**2).sum() / (len(ydata)- len(variables))
In (untested!) code, it might look like:
bestvals, covar, info, errmsg, ier = leastsq(objectivefunc, initvals,
args=args, full_output=1, **kws)
residual = objectivefunc(bestvals, *args)
reduced_chi_square = (residual**2).sum() / (len(data) - len(initvals))
covar = covar * reduced_chi_square
std_error = array([sqrt(covar[i, i]) for i in range(len(bestvals))])
Here, objectivefunc is expected to return the scaled misfit array,
that is "(data - model)/uncertainty".
Cheers,
--Matt Newville <newville at cars.uchicago.edu> 630-252-0431
More information about the SciPy-User
mailing list