How do I include errors for my data in the lmfit least squares miniimization, and what is this error for conf_interval2d function in lmfit?

Question

I am new to python, and trying to use the lmfit package to check my own calculations, however I am unsure (1) as to how to include the errors for data (sig) for the following test (and 2) of an error I get with conf_interval2d shown below):

    import numpy as np
    from lmfit import Parameters, Minimizer, conf_interval, conf_interval2d, minimize, printfuncs


    x=np.array([ 0.18,  0.26,  1.14,  0.63,  0.3 ,  0.22,  1.16,  0.62,  0.84,0.44,  1.24,  0.89,  1.2 ,  0.62,  0.86,  0.45,  1.17,  0.59, 0.85,  0.44])
    data=np.array([ 68.59,  71.83,  22.52,44.587,67.474 ,  55.765,  20.9,41.33783784,45.79 ,  47.88,   6.935,  34.15957447,44.175,  45.89230769,  57.29230769,  60.8,24.24335594,  34.09121287,  42.21504003,  26.61161674])
    sig=np.array([ 11.70309409,  11.70309409,  11.70309409,  11.70309409,11.70309409,  11.70309409,  11.70309409,  11.70309409,11.70309409,  11.70309409,  11.70309409,  11.70309409,11.70309409,  11.70309409,  11.70309409,  11.70309409,11.70309409,  11.70309409,  11.70309409,  11.70309409])

    def residual(pars, x, data=None):
        a=pars['a'].value
        b=pars['b'].value
        model = a + (b*x)
        if data is None:
            return model
        return model-data

    params=Parameters()
    params.add('a', value=70.0)
    params.add('b', value=40.0)

    mi=minimize(residual, params, args=(x, data))
    #mi=minimize(residual, params, args=(x,), kws={'data': data})#is this more correct?
    ci, trace = conf_interval(mi, trace=True)

This works fine until now, but as questioned above, how do I include the errors for data (sig_chla) so that I can calculated a weighted and reduced chi-square?

Part 2: FURTHERMORE, when I attempt to use the following so that I can plot the confidence intervals, xs, ys, grid = conf_interval2d(mi, 'a', 'b', 20, 20)

I get the following error:

* ValueError: failed to create intent(cache|hide)|optional array-- must have defined dimensions but got (0,)

or

Parameter 4 to routine DGESV was incorrect Mac OS BLAS parameter error in DGESV , parameter #0, (unavailable), is 0

Answer 1

You must weight your data by the errors yourself in the residual() function.

From the lmfit docs (not very easy to find, though):

Note that the calculation of chi-square and reduced chi-square assume that the returned residual function is scaled properly to the uncertainties in the data. For these statistics to be meaningful, the person writing the function to be minimized must scale them properly.

It is not so hard to do, though. From the Wikipedia entry on weighted least-square fitting :

If, however, the measurements are uncorrelated but have different uncertainties, a modified approach might be adopted. Aitken showed that when a weighted sum of squared residuals is minimized, is BLUE if each weight is equal to the reciprocal of the variance of the measurement.

However, lmfit takes in the residual, not the squared residual, so instead of just going

    # This is what you do with no errorbars -> equal weights.
    resids = model - data
    return resids

you should do something like this (with scipy imported as sp):

    # Do this to include errors as weights.
    resids = model - data
    weighted = sp.sqrt(resids ** 2 / sig ** 2)
    return weighted

This should give you the correctly weighted fit.