Python lmfit reduced chi-square too small after weighted fit

Question

I am running a fit in Python 2.7 with lmfit using some test data with the following code. I require a weighted fit with weights of 1/y (with the Leven-Marq. routine). I have defined weights and am using them here:

from __future__ import division
from numpy import array, var
from lmfit import Model
from lmfit.models import GaussianModel, LinearModel

import matplotlib.pyplot as plt
import seaborn as sns

xd = array([1267, 1268, 1269, 1270, 1271, 1272, 1273, 1274, 1275, 1276,
    1277, 1278, 1279, 1280, 1281, 1282, 1283, 1284, 1285, 1286, 1287, 1288,
     1289, 1290, 1291, 1292, 1293, 1294, 1295, 1296, 1297, 1298, 1299, 1300,
     1301, 1302, 1303, 1304, 1305, 1306, 1307, 1308, 1309, 1310, 1311, 1312,
     1313, 1314, 1315, 1316, 1317, 1318, 1319, 1320, 1321, 1322, 1323, 1324,
     1325, 1326, 1327, 1328, 1329, 1330, 1331, 1332, 1333, 1334])
yd = array([238, 262, 255, 271, 270, 281, 261, 278, 280, 254, 289, 285, 304, 314,
    329, 342, 379, 450, 449, 564, 613, 705, 769, 899, 987, 1043, 1183, 1295, 1298,
    1521, 1502, 1605, 1639, 1572, 1659, 1558, 1476, 1397, 1267, 1193, 1016, 951,
    835, 741, 678, 558, 502, 480, 442, 399, 331, 334, 308, 283, 296, 265, 264, 
    273, 258, 270, 262, 263, 239, 263, 251, 246, 246, 234])

mod = GaussianModel() + LinearModel()
pars  = mod.make_params(amplitude=25300, center=1299, sigma=7, slope=0, intercept=450)
result = mod.fit(yd, pars, method='leastsq', x=xd, weights=1./yd)
rsq = 1 - result.residual.var() / var(yd)
print(result.fit_report())
print rsq

plt.plot(xd, yd,         'bo', label='raw')
plt.plot(xd, result.init_fit, 'k--', label='Initial_Guess')
plt.plot(xd, result.best_fit, 'r-', label='Best')
plt.legend()
plt.show()

The output is:

[[Model]]
    (Model(gaussian) + Model(linear))
[[Fit Statistics]]
    # function evals   = 27
    # data points      = 68
    # variables        = 5
    chi-square         = 0.099
    reduced chi-square = 0.002
    Akaike info crit   = -434.115
    Bayesian info crit = -423.017
[[Variables]]
    sigma:       7.57360038 +/- 0.063715 (0.84%) (init= 7)
    center:      1299.41410 +/- 0.071046 (0.01%) (init= 1299)
    amplitude:   25369.3304 +/- 263.0961 (1.04%) (init= 25300)
    slope:      -0.15015228 +/- 0.071540 (47.65%) (init= 0)
    intercept:   452.838215 +/- 93.28860 (20.60%) (init= 450)
    fwhm:        17.8344656 +/- 0.150037 (0.84%)  == '2.3548200*sigma'
    height:      1336.33919 +/- 17.28192 (1.29%)  == '0.3989423*amplitude/max(1.e-15, sigma)'
.
.
.
.
0.999999993313

The last line (just above here, or immediately before plt.plot(xd, yd, 'bo', label='raw')) is the R^2 and the resulting fit is attached here..

The R^2 and the visual inspection of the output suggest this is a reasonable fit. I am expecting a reduced chi-squared of order 1.00 (source). However, the returned value for the reduced chi-squared value is several orders of magnitude smaller than 1.00.

Since the default is no weights in lmfit and I need a weighted fit, I have defined weights, but I think I need to be specifying them differently. My suspicion is this specification of weights might be causing the reduced chi-squared to be so small.

Is there a different way to specify weights, or some other parameter, such that the reduced chi-squared after the curve fit is close to or on the same magnitude as 1.00?

score 5 · Accepted Answer · answered Apr 15 '17 at 13:56

5

The weight in lmfit is a multiplicative factor for the residua to be minimized in the least-squares sense. That is, it replaces

residual = model - data

with

residual = (model - data) * weights

A common approach, and one that I think you might be intending, is to say that the weights should be 1.0/variance_in_data, as that is what usually means to get to reduced chi-square around 1 for a good fit, as the excellent writeup you link to discusses.

As discussed there, the problem is determining the variance in the data. For many cases, such as when the signal is dominated by counting statistics, the variance in data can be estimated as sqrt(data). This ignores many sources of noise, but is often a good starting point. As it happens, I believe using

result = model.fit(..., weights=np.sqrt(1.0/yd))

will lead to a reduced chi-square of around 0.8 for your case. I think that is probably what you want.

Also, to clarify a related point: The writeup you link discusses scaling the uncertainties in the fitted parameters when reduced chi-square is far from 1. Lmfit does the scaling described there by default (the scale_covar option can turn this off), so that changing the scale of the weights won't change the scale of the uncertainties in the parameters sigma, center, etc. The values for the uncertainties (and, best-fit values) will change some because the change in weighting changes the emphasis for each data point, but the best-fit values won't change much, and the estimated uncertainties should stay the same order of magnitude even if your estimate of the variance in the data (and so reduced chi-square) is off by a few orders of magnitude.

That is, changing your script to use weights=1.0/np.sqrt(yd) will bring reduced chi-square much closer to 1, but it will not change the uncertainties in the fitted variables very much.

answered Apr 15 '17 at 13:56

M Newville

7,486
2
16
29

Okay three questions: I indeed got the reduced chi-square to be `0.790` with your specification for weights. I tried a manual calculation of reduced chi-square with `np.sum(((yd - result.best_fit)**2)/result.best_fit)/(68-5)` and got a slightly different value of `0.78065`. I took 68 as the number of points and 5 as the number of fitting parameters (i.e. number of constraints in the problem) so the number of degrees of freedom is 63. The difference is almost negligible....does `lmfit` use a slightly different approach to estimating reduced chi-sq.? – edesz Apr 15 '17 at 17:51
Thanks for the great answer! Another follow-up question: I didn't think of this until your comment about scaling uncertainties but: when I use (a) `scale_covar=True`, the parameter error (for example) in amplitude is closer to the 67.4% confidence level reported from `result.ci_report()` versus (b) `scale_covar=False`, the parameter error in amplitude is considerably larger than the 67.4% confidence level reported from `result.ci_report()`. If 1*`sigma` par. uncertainties are wanted, should one only use the `.ci_report()` output? Is there a way to get `.fit_report()` uncertainties at 1`sigma`? – edesz Apr 15 '17 at 18:33
Regarding the weights - yes, using `sqrt(yd)` gives the red-chi sq. values close to 1.0. Thanks again for the great explanation. Okay, my 3rd and the main question I have about your answer: you have referred to variance as being considered `sqrt(yd)`. My understanding of weights([wiki link](https://en.wikipedia.org/wiki/Least_squares#Weighted_least_squares)) is they should be the reciprocal of variance (i.e. weights = 1 / `sigma`^2) where variance = `sigma`^2. You have also used reciprocal of variance, but is `sqrt(yd)` referring to standard deviation (`sigma`, i.e. square root of variance)? – edesz Apr 15 '17 at 18:57
1

The reduced chi-square should indeed be `chi-square/63` for your case. But `chi-square` is the sum of squares of the residual array which includes the weighting factor. For Question 2: the uncertainties reported are meant to be the 1-sigma uncertainties, and are automatically calculated based on the fit process. Compared to the explicit, brute-force `conf_intervals` approach, the automatic errors should be viewed as approximate, but for your fit are probably really close. And a meta-coment: if the answer helps, the SO convention is to accept and/or upvote it. – M Newville Apr 16 '17 at 01:47
1

The third question is essentially "what is the variance in a signal", which depends greatly on the nature of the signal. If the signal is "counts" of something, then "shot noise" or "counting statistics" often dominates, and the variance goes as sqrt(signal). For many physical processes, this is a decent starting point -- it may ignore other contributions to the noise ("systematic errors"), but these are usually harder to assess without detailed understanding of the underlying nature of the signal and its measurements. – M Newville Apr 16 '17 at 01:52
Thanks for all the info. I'm still thinking about `scale_covar` since your post here and I *will* ask a separate question about that in a few days. As far as this question goes about weighting, and follow-ups, everything seems good. No further questions. Thanks again....great info. – edesz Apr 25 '17 at 20:41
@MNewville sorry if I revive this old thread, I'm comparing some fitting results between `lmfit` and `ROOT`, because I'd like to use Python but everyone trusts `ROOT` to do the right thing out of the box, especially when it comes to fit errors estimates. I got to a point where I get similar results, except `ROOT` scales fit errors with data errors so I need to set `scale_covar=False` in `lmfit` to get comparable estimates. Intuitively that seems correct, I'd expect bigger errors if data is more noisy but there must also be a good reason `lmfit` scales covariance by default. Who should I trust? – filippo Dec 18 '18 at 18:04
@MNewville I see the issue has been discussed plenty of times in the mailing-list, will do some more digging there. Thanks for the great module! – filippo Dec 18 '18 at 18:12
@filippo Using `scale_covar=True` is basically admitting that you don't know the scale of the uncertainties in the data very well, but do assume that the current fit is "good" (otherwise, you probably aren't going to use the estimated uncertainties in the parameters anyway). That is, it reports uncertainties as if `epsilon` were set so that reduced chi-square = 1. In many real measurements there is not a very good estimate of data uncertainties, so this is useful. But if you have a good estimate of data uncertainties, by all means use those and set `scale_covar=False`. – M Newville Dec 18 '18 at 20:21
@MNewville thank you much clearer now. It's mostly histograms from spectroscopy with `sqrt(N)` uncertainties. The thing is parameters with estimated uncertainties are being used in subsequent weighted fits and somehow with `scale_covar=True` this leads to somewhat optimistic uncertainties, at least compared to the same procedure done in Root. Still trying to better understand it all, as Root documentation can sometimes be a little obscure. – filippo Dec 18 '18 at 21:02
@filippo it's a little hard to have a detailed discussion in a comments section. Posting a question on SO or using the lmfit mailing list would allow you to give a more complete idea of what you're trying to do and what you find for differences in the results from root vs lmfit. – M Newville Dec 18 '18 at 21:08
@MNewville right, I'll try to come up with a minimal test case and ask in the mailing list. Thanks again for now. – filippo Dec 18 '18 at 21:13
In case other [SO](https://stackoverflow.com/) users come here and find that a similar discussion would also be useful to them, here is the direct link to the `lmfit` mailing list: [link](https://groups.google.com/forum/#!forum/lmfit-py) – edesz Dec 19 '18 at 02:42
@filippo I enjoyed your discussion and I wonder if you found your answer. – Delosari Oct 14 '20 at 12:11
@Delosari nothing really conclusive, my experience is that most people trust the library they're using for fits without further question. In physics you generally have a good estimate of experimental uncertainties and you expect big fit errors if you have big experimental error bars. Physics oriented libs, like ROOT, behave this way. Statistics oriented libs like most in the python world treat uncertainties as relative weights and scale covariance matrix to reflect this. Problems arise when you use a mix of both, and you need to tune your fit libraries behave the same way. – filippo Oct 14 '20 at 12:42
@Delosari my advice as a rule of thumb would be: find what the most accepted software in your field does, and tune your python libraries to do the same :-D – filippo Oct 14 '20 at 12:44
@Delosari note that this goes beyond lmfit itself, most python libraries that do some kind of fitting have a parameter to configure this particular aspect and they usually default to relative weights and covariance scaling. – filippo Oct 14 '20 at 12:49
@filippo :D Thank you that is a wise advice – Delosari Oct 14 '20 at 13:05

Python lmfit reduced chi-square too small after weighted fit

1 Answers1

Linked