2017-04-15 14 views
1

다음 코드로 일부 테스트 데이터를 사용하여 Python 2.7에서 lmfit으로 적합하게 실행 중입니다. 1/y (Leven-Marq 루틴과 함께)의 가중치가 맞는 것을 요구합니다. 나는 가중치를 정의하고 여기에 사용하고 있습니다 :Python lmfit 가중치 적용 후 카이 제곱이 너무 작음

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 및 결과 적합 여기에 첨부되어 있습니다. enter image description here.

R^2와 출력의 육안 검사는 이것이 합리적인 적합이라고 제안합니다. 주문 1.00 (source)의 축소 된 카이 제곱을 기대합니다. 그러나 축소 된 카이 제곱 값에 대한 반환 값은 1.00보다 몇 배 작은 값입니다.

lmfit의 기본값은 no weights이며 가중치가 필요합니다. 가중치를 정의했지만,이를 다르게 지정해야한다고 생각합니다. 저의 의혹은이 가중치의 지정이 축소 된 카이 제곱을 너무 작게 만들 수 있습니다.

커브 피팅 후 축소 된 카이 제곱이 1.00에 근접하거나 같은 크기가되도록 가중치 또는 다른 매개 변수를 지정하는 다른 방법이 있습니까?

답변

2

lmfit의 가중치는 최소 자승 감에서 최소화되는 잔류 물에 대한 곱셈 적 요소입니다. 즉 같은 즉, 난 당신이 가중치는 1.0/variance_in_data해야한다는 말을하는 것입니다, 의도 것 같아요

residual = (model - data) * weights 

일반적인 방법, 그리고 하나

residual = model - data 

을 대체 무엇인지 일반적으로 훌륭한 적합성을 위해 축소 된 카이 제곱에 도달하는 것을 의미합니다.

거기에서 설명한대로 문제는 데이터의 분산을 결정하는 것입니다. 신호가 계수 통계에 의해 지배되는 경우와 같이 많은 경우에 데이터의 분산은 sqrt(data)으로 추정 할 수 있습니다. 이것은 많은 소음원을 무시하지만 종종 좋은 출발점입니다. 결과적으로 나는 사용을 믿습니다

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

귀하의 경우 약 0.8의 카이 제곱으로 이어질 것입니다. 아마 그게 네가 원하는 것 같아. 또한

는 관련 점을 명확히 : 카이 제곱 감소 할 때 링크를 작성자가 장착 매개 변수의 불확실성을 확장에 대해 설명 Lmfit 기본적으로이 기술 스케일링합니다 ( scale_covar 옵션이 해제 할 수 있습니다 않습니다 1. 거리가 멀다) 따라서 가중치의 변경은 매개 변수 sigma, center 등의 불확실성의 크기를 변경하지 않습니다. 불확실성 (및 최적치)에 대한 값은 가중치 변경의 변경으로 인해 일부 변경됩니다 각 데이터 포인트에 대한 강조점이지만 가장 적합한 값은 많이 변경되지 않으며 데이터의 분산 추정치 (예 : reduced chi-square)가 몇 개씩 떨어져도 추정 된 불확실성은 동일한 크기로 유지되어야합니다 규모의 순서.

즉, weights=1.0/np.sqrt(yd)을 사용하도록 스크립트를 변경하면 카이 제곱이 1에 훨씬 가까워 지지만 적합 변수의 불확실성이 크게 변경되지는 않습니다.

+0

좋아요 세 가지 질문 : 실제로 무게에 대한 사양이 '0.790'으로 줄어 들었습니다. 나는'np.sum ((yd - result.best_fit) ** 2) /result.best_fit)/(68-5)'를 사용하여 축소 된 카이 제곱의 수동 계산을 시도했으며'0.78065'의 값이 약간 다릅니다 . 필자는 68 점을 피팅 매개 변수의 수 (즉, 문제의 제약 수)로 취하여 자유도의 수는 63입니다. 차이는 거의 무시할 수 있습니다 .... lmfit은 감소 된 카이 스퀘어를 추정하기위한 약간 다른 접근법? –

+0

대단한 답변을 보내 주셔서 감사합니다! 또 다른 후속 질문 : 불확실성을 조정하는 것에 대해 의견을 말할 때까지는 생각하지 않았습니다 : (a)'scale_covar = True'를 사용하면 진폭의 매개 변수 오류 (예 :)가 67.4 %의 신뢰 수준에 더 가깝습니다 'result.ci_report()'에서 (b)'scale_covar = False'로보고되면, 진폭의 매개 변수 오류는'result.ci_report()'에서보고 된 67.4 % 신뢰 수준보다 상당히 큽니다. 1 *'시그마'파. '.ci_report()'출력 만 사용해야한다면 불확실성이 필요합니까? 1' 시그마에서 .fit_report() 불확실성을 얻을 수있는 방법이 있습니까? –

+0

가중치 - 예,'sqrt (yd) '를 사용하면 적색 - 카이 평방 값이 1.0에 가깝습니다. 훌륭한 설명에 다시 한 번 감사드립니다. 좋아요, 제 3 그리고 당신의 대답에 대한 주요 질문 : 분산을 sqrt (yd)로 간주했습니다. 가중치에 대한 나의 이해 ([wiki link] (https://en.wikipedia.org/wiki/Least_squares#Weighted_least_squares))는 분산의 역수 여야합니다 (예 : 가중치 = 1 /'시그마^2). 여기서 variance = '시그마'^ 2. 또한 분산의 역수를 사용했지만 표준 편차 ('시그마 ', 즉 제곱근의 분산)를 나타내는'sqrt (yd)'입니까? –