Python lmfit加权拟合后减小的卡方值过小

Question

Python lmfit加权拟合后减小的卡方值过小

python-2.7curve-fittingleast-squareslmfit

4

我正在使用以下代码在Python 2.7中运行lmfit拟合一些测试数据。我需要进行加权拟合，并使用1/y作为权重（使用Leven-Marq.算法）。我已经定义了权重并在此处使用：

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()

输出为：

[[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

在这里（或者紧接着 plt.plot(xd, yd, 'bo', label='raw') 之前）的最后一行是 R^2，生成的拟合结果附在这里。

R^2 和输出的可视检查表明这是一个合理的拟合。我预计减小的卡方应该为 1.00 的量级（来源）。然而，返回的减小的卡方值比 1.00 小了几个数量级。

由于在 lmfit 中默认是不使用权重，而我需要进行加权拟合，所以我定义了权重，但我认为我需要以不同的方式指定它们。我怀疑这种权重的规范可能导致减小的卡方值太小。

有没有其他方式来指定权重或其他参数，使得曲线拟合后的降低卡方接近或与1.00相同的数量级？

- edesz

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- M Newville · Accepted Answer

lmfit中的权重是使残差在最小二乘意义下进行缩放的乘法因子。也就是说，它取代了

residual = model - data

使用

residual = (model - data) * weights

一种常见的方法是将权重设置为1.0/数据方差，这通常可以实现较好的拟合效果，使得缩小后的卡方值接近于1。正如您提供的链接所述，这是一个很好的解释。

但是，确定数据的方差是一个问题。对于许多情况，例如在信号被计数统计占主导地位时，数据的方差可以估计为sqrt(data)。虽然这忽略了许多噪声源，但通常可以作为一个良好的起点。恰好，我相信使用

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

对于您的情况，这将导致减少卡方约为0.8。我想这可能是您想要的。

另外，为了澄清一个相关点：您链接的写作讨论了当减少的卡方远离1时缩放拟合参数的不确定性。Lmfit默认执行了描述在那里缩放（scale_covar选项可以关闭此功能），因此改变权重的比例不会改变参数sigma、center等的不确定度的比例。不确定度的值（以及最佳拟合值）会有所改变，因为加权的变化改变了每个数据点的重点，但最佳拟合值不会有太大变化，即使您对数据的方差估计（因此reduced chi-square）差了几个数量级，估计的不确定度也应该保持相同的数量级。

也就是说，将您的脚本更改为使用weights=1.0/np.sqrt(yd)将使减少的卡方更接近1，但它不会对拟合变量的不确定度产生太大影响。