[SciPy-User] optimize.leastsq error estimation

[Matt Newville]

Hi Joe,

On Thu, Mar 21, 2013 at 3:10 PM, Joe Philip Ninan <indiajoe at 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