4

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

Nicole Goebel
  • 537
  • 1
  • 9
  • 17

1 Answers1

6

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.

Thriveth
  • 377
  • 1
  • 6
  • 18
  • What does the error mean here? is it a variance? a standard deviation?? – Mehdi Oct 16 '15 at 09:05
  • Sorry, "error" in this case means standard deviation. – Thriveth Oct 16 '15 at 13:02
  • I don't think this is right as the docs explicitly says that the residual function must return the residuals only. – rhody Aug 07 '16 at 22:00
  • @rhody The residuals are returned for minimization purposes only. You are free to return any array for minimization - and the above described procedure gives the statistically correct minimization vector for obtaining the weighted fit. – Thriveth Aug 08 '16 at 14:47
  • You're right, I found another sentence buried in the docs that says we can weight the residuals which I've done and it works. However I am not sure if one should compute the full chi-squared in tehe residuals function. The docs aren't clear here but I am assuming that after the call to residuals the lmfit itself computes the chi-squared? – rhody Aug 09 '16 at 16:37
  • That is a good point, and that is exactly why I let the code in the example return the "weighted-chi-non-squared". – Thriveth Aug 09 '16 at 22:26
  • Can you explain again, why not `weighted = resids / sig` but `weighted = sp.sqrt(resids ** 2 / sig ** 2)` ?? After almost 10 years :) – KaO Dec 01 '22 at 20:00