Python lmfit在加權合適後縮小了卡方太小

Question

我正在使用一些測試數據並使用以下代碼在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的訂單減少卡方（source）。但是，降低的卡方值的返回值比1.00小几個數量級。

由於默認值是no weightslmfit我需要一個加權擬合，我已經定義了權重，但我認爲我需要以不同的方式指定它們。我的懷疑是這種重量的規格可能會導致減小的卡方非常小。

是否有不同的方式來指定權重或其他參數，使得曲線擬合後的減少的卡方接近或等於1.00的相同量級？

Answer 1

lmfit中的權重是殘差在最小二乘意義上被最小化的乘法因子。也就是說，它取代

residual = model - data 

與

residual = (model - data) * weights 

一種常見的做法，一個是我想你可能會打算，是說的權重應爲1.0/variance_in_data，因爲這是通常意味着要在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，但它不會很大程度上改變擬合變量的不確定性。