python - LMFIT 曲线拟合模块有时不显示错误
问题描述
我有各种试图拟合函数的频谱图。LMFIT 有一个我正在使用的复合模型功能，我的模型本质上是在恒定背景下的 Voigt 或高斯峰的总和。峰值中心的初始猜测是使用 scipy peak finder 功能找到的。
拟合实际上大部分都非常好，即使对于具有许多峰值的数据集也是如此。我的问题是有时（对于一些具有更多峰值的较大数据集）拟合报告上的错误不会显示，但有时会显示。
通过阅读有关类似主题的其他问题，我似乎已经收集到这可能是由于未使用参数，或者被推到我设置的边界。鉴于此，我尝试删除所有边界，并为每个峰值设置初始条件。它仍然没有显示错误。
我意识到将 15 条左右重叠的高斯曲线拟合到频谱图可能需要很多，但由于肉眼看来拟合效果很好，我认为它必须在某种程度上起作用。
本质上，我想知道如何在拟合报告中为具有许多峰值的数据集获取错误，就像具有 2 个峰值的数据集一样。
以下是两个示例拟合报告：
两个峰值，显示错误：
[[Model]]
((Model(constant) + Model(voigt, prefix='g0_')) + Model(voigt, prefix='g1_'))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 186
# data points = 564
# variables = 9
chi-square = 1342.58147
reduced chi-square = 2.41906571
Akaike info crit = 507.154526
Bayesian info crit = 546.170014
[[Variables]]
g0_sigma: 7.58224056 +/- 0.07787554 (1.03%) (init = 7)
g0_center: 657.390036 +/- 0.02496585 (0.00%) (init = 656.4029)
g0_amplitude: 1549654.58 +/- 5871.81065 (0.38%) (init = 57643.33)
g0_gamma: 3.70825451 +/- 0.10708503 (2.89%) (init = 0.7)
g0_fwhm: 27.3059988 +/- 0.28045426 (1.03%) == '3.6013100*g0_sigma'
g0_height: 57414.5545 +/- 137.261488 (0.24%) == 'g0_amplitude*wofz((1j*g0_gamma)/(g0_sigma*sqrt(2))).real/(g0_sigma*sqrt(2*pi))'
g1_sigma: 7.66546221 +/- 0.55266628 (7.21%) (init = 7)
g1_center: 803.744461 +/- 0.16381855 (0.02%) (init = 803.2903)
g1_amplitude: 315684.011 +/- 6421.21329 (2.03%) (init = 10676.22)
g1_gamma: 6.21905616 +/- 0.67178258 (10.80%) (init = 0.7)
g1_fwhm: 27.6057057 +/- 1.99032260 (7.21%) == '3.6013100*g1_sigma'
g1_height: 9525.49591 +/- 130.820663 (1.37%) == 'g1_amplitude*wofz((1j*g1_gamma)/(g1_sigma*sqrt(2))).real/(g1_sigma*sqrt(2*pi))'
c: 1096.67100 +/- 6.42230443 (0.59%) (init = 0)
[[Correlations]] (unreported correlations are < 0.250)
C(g0_sigma, g0_gamma) = -0.929
C(g1_sigma, g1_gamma) = -0.926
C(g0_amplitude, g0_gamma) = 0.828
C(g1_amplitude, g1_gamma) = 0.821
C(g0_sigma, g0_amplitude) = -0.659
C(g1_sigma, g1_amplitude) = -0.650
然后没有错误（忽略减少的卡方等）：
[[Model]]
((((((((((((((((((Model(constant) + Model(voigt, prefix='g0_')) + Model(voigt, prefix='g1_')) + Model(voigt, prefix='g2_')) + Model(voigt, prefix='g3_')) + Model(voigt, prefix='g4_')) + Model(voigt, prefix='g5_')) + Model(voigt, prefix='g6_')) + Model(voigt, prefix='g7_')) + Model(voigt, prefix='g8_')) + Model(voigt, prefix='g9_')) + Model(voigt, prefix='g10_')) + Model(voigt, prefix='g11_')) + Model(voigt, prefix='g12_')) + Model(voigt, prefix='g13_')) + Model(voigt, prefix='g14_')) + Model(voigt, prefix='g15_')) + Model(voigt, prefix='g16_')) + Model(voigt, prefix='g17_'))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 101758
# data points = 564
# variables = 73
chi-square = 13631.1513
reduced chi-square = 27.7620190
Akaike info crit = 1942.37313
Bayesian info crit = 2258.83209
[[Variables]]
g0_sigma: 192.689563 (init = 7)
g0_center: 422.997773 (init = 389.2612)
g0_amplitude: 1068.20554 (init = 2820.275)
g0_gamma: -618.292505 (init = 0.7)
g0_fwhm: 693.934849 == '3.6013100*g0_sigma'
g0_height: 760.695408 == 'g0_amplitude*wofz((1j*g0_gamma)/(g0_sigma*sqrt(2))).real/(g0_sigma*sqrt(2*pi))'
g1_sigma: 17.9116349 (init = 7)
g1_center: 431.473501 (init = 431.5196)
g1_amplitude: 4525.79900 (init = 3929.55)
g1_gamma: -36.4029462 (init = 0.7)
g1_fwhm: 64.5053499 == '3.6013100*g1_sigma'
g1_height: 1556.62320 == 'g1_amplitude*wofz((1j*g1_gamma)/(g1_sigma*sqrt(2))).real/(g1_sigma*sqrt(2*pi))'
g2_sigma: 4.86138247 (init = 7)
g2_center: 805.214348 (init = 803.2903)
g2_amplitude: 668696.572 (init = 33620)
g2_gamma: 3.52118850 (init = 0.7)
g2_fwhm: 17.5073453 == '3.6013100*g2_sigma'
g2_height: 33446.8793 == 'g2_amplitude*wofz((1j*g2_gamma)/(g2_sigma*sqrt(2))).real/(g2_sigma*sqrt(2*pi))'
g3_sigma: 5.13832814 (init = 7)
g3_center: 1032.12566 (init = 1035.003)
g3_amplitude: 595227.235 (init = 17401.5)
g3_gamma: 8.79672669 (init = 0.7)
g3_fwhm: 18.5047125 == '3.6013100*g3_sigma'
g3_height: 17386.9619 == 'g3_amplitude*wofz((1j*g3_gamma)/(g3_sigma*sqrt(2))).real/(g3_sigma*sqrt(2*pi))'
g4_sigma: 5.06870799 (init = 7)
g4_center: 1160.54035 (init = 1160.308)
g4_amplitude: 74021.3519 (init = 4387.175)
g4_gamma: 5.14599307 (init = 0.7)
g4_fwhm: 18.2539888 == '3.6013100*g4_sigma'
g4_height: 3023.66817 == 'g4_amplitude*wofz((1j*g4_gamma)/(g4_sigma*sqrt(2))).real/(g4_sigma*sqrt(2*pi))'
g5_sigma: 5.71945320 (init = 7)
g5_center: 1270.96477 (init = 1268.509)
g5_amplitude: 458428.440 (init = 15159)
g5_gamma: 7.80744109 (init = 0.7)
g5_fwhm: 20.5975240 == '3.6013100*g5_sigma'
g5_height: 13982.1731 == 'g5_amplitude*wofz((1j*g5_gamma)/(g5_sigma*sqrt(2))).real/(g5_sigma*sqrt(2*pi))'
g6_sigma: 2.6981e-09 (init = 7)
g6_center: 1448.28921 (init = 1352.596)
g6_amplitude: 572109.865 (init = 4058.475)
g6_gamma: 13.5568440 (init = 0.7)
g6_fwhm: 9.7166e-09 == '3.6013100*g6_sigma'
g6_height: 13432.9366 == 'g6_amplitude*wofz((1j*g6_gamma)/(g6_sigma*sqrt(2))).real/(g6_sigma*sqrt(2*pi))'
g7_sigma: 12.5995161 (init = 7)
g7_center: 1351.19388 (init = 1450.939)
g7_amplitude: 64633.5688 (init = 13943)
g7_gamma: -1.76986853 (init = 0.7)
g7_fwhm: 45.3747632 == '3.6013100*g7_sigma'
g7_height: 2297.69041 == 'g7_amplitude*wofz((1j*g7_gamma)/(g7_sigma*sqrt(2))).real/(g7_sigma*sqrt(2*pi))'
g8_sigma: 697.890539 (init = 7)
g8_center: 2422.51062 (init = 1855.019)
g8_amplitude: 210.220409 (init = 2548.15)
g8_gamma: -3083.15825 (init = 0.7)
g8_fwhm: 2513.32018 == '3.6013100*g8_sigma'
g8_height: 4158.40426 == 'g8_amplitude*wofz((1j*g8_gamma)/(g8_sigma*sqrt(2))).real/(g8_sigma*sqrt(2*pi))'
g9_sigma: 83.3565885 (init = 7)
g9_center: 2867.11059 (init = 1926.474)
g9_amplitude: 6021213.64 (init = 2513.125)
g9_gamma: -231.269811 (init = 0.7)
g9_fwhm: 300.192916 == '3.6013100*g9_sigma'
g9_height: 2697769.12 == 'g9_amplitude*wofz((1j*g9_gamma)/(g9_sigma*sqrt(2))).real/(g9_sigma*sqrt(2*pi))'
g10_sigma: 105.224438 (init = 7)
g10_center: 2757.10276 (init = 1940.692)
g10_amplitude: -13348456.2 (init = 2628.65)
g10_gamma: -300.405844 (init = 0.7)
g10_fwhm: 378.945820 == '3.6013100*g10_sigma'
g10_height: -5945307.23 == 'g10_amplitude*wofz((1j*g10_gamma)/(g10_sigma*sqrt(2))).real/(g10_sigma*sqrt(2*pi))'
g11_sigma: 100.413181 (init = 7)
g11_center: 2741.65736 (init = 2226.984)
g11_amplitude: 10844230.4 (init = 2359.375)
g11_gamma: -296.216143 (init = 0.7)
g11_fwhm: 361.618992 == '3.6013100*g11_sigma'
g11_height: 6673387.09 == 'g11_amplitude*wofz((1j*g11_gamma)/(g11_sigma*sqrt(2))).real/(g11_sigma*sqrt(2*pi))'
g12_sigma: 366.846051 (init = 7)
g12_center: 1940.87402 (init = 2450.446)
g12_amplitude: 256322.806 (init = 2311.475)
g12_gamma: -753.396963 (init = 0.7)
g12_fwhm: 1321.12635 == '3.6013100*g12_sigma'
g12_height: 4501.31993 == 'g12_amplitude*wofz((1j*g12_gamma)/(g12_sigma*sqrt(2))).real/(g12_sigma*sqrt(2*pi))'
g13_sigma: 103.888422 (init = 7)
g13_center: 3019.80754 (init = 2667.96)
g13_amplitude: 1616371.56 (init = 4225.95)
g13_gamma: -375.306839 (init = 0.7)
g13_fwhm: 374.134413 == '3.6013100*g13_sigma'
g13_height: 8468440.28 == 'g13_amplitude*wofz((1j*g13_gamma)/(g13_sigma*sqrt(2))).real/(g13_sigma*sqrt(2*pi))'
g14_sigma: 147.814985 (init = 7)
g14_center: 2947.03711 (init = 2700.413)
g14_amplitude: 167.834745 (init = 3302.95)
g14_gamma: -836.092487 (init = 0.7)
g14_fwhm: 532.327584 == '3.6013100*g14_sigma'
g14_height: 8027180.83 == 'g14_amplitude*wofz((1j*g14_gamma)/(g14_sigma*sqrt(2))).real/(g14_sigma*sqrt(2*pi))'
g15_sigma: 145.620972 (init = 7)
g15_center: 2989.06096 (init = 2803.39)
g15_amplitude: -886458.453 (init = 2653.55)
g15_gamma: -516.771454 (init = 0.7)
g15_fwhm: 524.426264 == '3.6013100*g15_sigma'
g15_height: -2636036.61 == 'g15_amplitude*wofz((1j*g15_gamma)/(g15_sigma*sqrt(2))).real/(g15_sigma*sqrt(2*pi))'
g16_sigma: 78.0169344 (init = 7)
g16_center: 3135.45756 (init = 2860.738)
g16_amplitude: -1388405.82 (init = 38420.5)
g16_gamma: -217.130421 (init = 0.7)
g16_fwhm: 280.963166 == '3.6013100*g16_sigma'
g16_height: -680872.187 == 'g16_amplitude*wofz((1j*g16_gamma)/(g16_sigma*sqrt(2))).real/(g16_sigma*sqrt(2*pi))'
g17_sigma: 81.7099268 (init = 7)
g17_center: 3322.62834 (init = 2936.564)
g17_amplitude: 333478.937 (init = 43586)
g17_gamma: -216.876994 (init = 0.7)
g17_fwhm: 294.262776 == '3.6013100*g17_sigma'
g17_height: 109848.342 == 'g17_amplitude*wofz((1j*g17_gamma)/(g17_sigma*sqrt(2))).real/(g17_sigma*sqrt(2*pi))'
c: 1502.83014 (init = 0)
这里还有一个精简版的代码（PseudoPseudoCode？），它没有运行，但显示了我正在尝试做的事情。我运行的实际代码所以它不是代码的真正问题，可以这么说，但也许它会说明我正在尝试做的事情。我敢肯定它已经过优化，所以请随意在你的回答中嘲笑我。
#read in data
#decide on stn ratio prominence and other values, the model type to use etc
#find estimate for peak centre, using this scipy function
peaks = scipy.signal.find_peaks(y,height = error*StN, width = width,prominence=PStN)
.
.
.
.
.
#obviously cut out some of the other code but this loops to make a gaussian or a voigt or whatever for each peak found above
model_array = {}
pars = lm.Parameters()
while n < peakamt:
model_array[n] = modeltype(prefix = "g"+str(n)+"_")
pars.update(model_array[n].make_params())
pars['g'+str(n)+'_center'].set(value=xpeaks[n])
pars['g'+str(n)+'_sigma'].set(value=7)
pars['g'+str(n)+'_amplitude'].set(value=ypeaks[n])
if modeltype == VoigtModel:
pars["g"+str(n)+'_gamma'].set(value=0.7, vary=True, expr = '')
#print(n)
n = n + 1
#adds the constant model
const = ConstantModel()
pars.update(const.make_params())
#makes the entire composite model, adding each peak model to the constant model
model = const
m = 0
while m < n:
model = model + model_array[m]
time.sleep(0.1)
print(".")
m = m+1
#here is the fitting
out = model.fit(y, pars, x=x, weights = weight)
final_fit = np.array(out.best_fit)
residuals = final_fit - y
#creating main figure
fig = plt.subplot(2,1,1)
#fig.plot(w,Lorentz)
fig.plot(x,final_fit)
fig.plot(x, y, 'ro', markersize = 2)
fig.errorbar(x,y,yerr=e, linestyle='none')
fig.set_title(str(ElementName) + " Deg " + str(Degrees))
fig.set_ylabel('Intensity (counts)')
fig.set_xlabel('Wavenumber (cm^-1)')
#creating residual figure
res = plt.subplot(2,1,2)
res.plot(x,residuals, 'r+')
res.errorbar(x,residuals,yerr=e,linestyle='none')
res.set_title("Residuals")
#Formatting the figures
plt.tight_layout()
#Output
print(out.fit_report(min_correl=0.25))
print('\n')
plt.show()
非常感谢您阅读本文，希望有人可以帮助我。如果这对阅读是一种折磨，我深表歉意，我对这个网站或一般编码没有太多经验。我希望它不会太糟糕。
解决方案
最新版本lmfit的报告中应该包含一些关于为什么它无法估计不确定性的信息。通常，如果某些参数偏离得太远以至于其贡献对模型或拟合没有影响，或者如果某些参数卡在边界上，就会发生这种情况。报告应该给出一些暗示。
在您的情况下，看起来gamma参数都变得非常负面并且完全疯狂。对于 Voigt 函数来说，gamma < 0并不完全是“荒谬的”，但它确实使该函数不完全是一个“峰值”，更像是“尖锐洛伦兹减去更宽高斯”的混合。您可能想为 设置一个下限gamma，也许是 0。
还有其他可能性，例如限制gamma为某些因子时间sigma并对所有峰值使用相同的因子，例如
pars = lm.Parameters()
pars.add('gamma_scale', value=0.7, min=0, max=100, vary=True)
while n < peakamt:
pref = 'g%d_' % n
model_array[n] = modeltype(prefix =pref)
pars.update(model_array[n].make_params())
pars['%s_center' % pref].set(value=xpeaks[n])
pars['%s_sigma' % pref].set(value=7)
pars['%s_amplitude' % pref].set(value=ypeaks[n])
if modeltype == VoigtModel:
pars['%s_gamma' % pref].set(expr='gamma_scale*%s_sigma' % pref)
n = n + 1
这仍将允许 Voigt 峰与 有一些可变性gamma，但将它们限制为都具有相同数量的gamma-ness。这是否是您真正想要做的可能取决于您的数据的性质。
推荐阅读
- javascript - 如何将值分配给 React.js 中另一个数组中的关联数组
- java - 如何突出显示选定的 RecyclerView 项目并在 Activity 中使用其值
- branch - 使用“if”语句进行分支
- azure - 从 Azure Function 访问 CosmosDB（无输入绑定）
- python - 堆栈是否在我的 python 程序中做一些有趣的事情
- javascript - RegExp 接受正数和负数
- python - Trying to convert Left Join of SQL into Django query-set?
- github - Initialized repo on the wrong folder
- swift - UIScrollView with vertical stackview scrolls horizontally but not vertically
- javascript - 使用错误的导航链接突出显示的 Scrollspy 活动类