LMFIT 曲线拟合模块有时不显示错误
LMFIT Curve Fitting Module Sometimes doesn't show Errors
我有各种各样的频谱图,我正试图将它们拟合到一个函数中。 LMFIT 具有我正在使用的复合模型功能,我的模型本质上是恒定背景下 Voigt 或高斯峰的总和。使用 scipy 峰值查找器功能找到峰值中心的初始猜测。
拟合实际上大部分都很好,即使对于具有许多峰的数据集也是如此。我的问题是有时(对于一些具有更多峰的较大数据集)拟合报告中的错误不会显示,但有时会显示。
通过阅读有关类似主题的其他问题,我似乎了解到这可能是由于未使用参数或被推到我设置的边界所致。鉴于此,我尝试移除所有边界,并为每个峰设置初始条件。它仍然没有显示错误。
我意识到将 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?),它没有 运行 但显示了我正在尝试做的事情。我有 运行s 的实际代码,所以这并不是代码的问题,可以这么说,但也许它会说明我正在尝试做的事情。我确定它的优化非常糟糕,所以请随意在你的回答中嘲笑我。
#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
不完全是 "absurd",但它确实使函数不完全是 "peak",更像是 "sharp Lorentzian minus broader Gaussian" 的混合。您可能想在 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。这是否是您真正想要做的可能取决于数据的性质。
我有各种各样的频谱图,我正试图将它们拟合到一个函数中。 LMFIT 具有我正在使用的复合模型功能,我的模型本质上是恒定背景下 Voigt 或高斯峰的总和。使用 scipy 峰值查找器功能找到峰值中心的初始猜测。
拟合实际上大部分都很好,即使对于具有许多峰的数据集也是如此。我的问题是有时(对于一些具有更多峰的较大数据集)拟合报告中的错误不会显示,但有时会显示。
通过阅读有关类似主题的其他问题,我似乎了解到这可能是由于未使用参数或被推到我设置的边界所致。鉴于此,我尝试移除所有边界,并为每个峰设置初始条件。它仍然没有显示错误。
我意识到将 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?),它没有 运行 但显示了我正在尝试做的事情。我有 运行s 的实际代码,所以这并不是代码的问题,可以这么说,但也许它会说明我正在尝试做的事情。我确定它的优化非常糟糕,所以请随意在你的回答中嘲笑我。
#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
不完全是 "absurd",但它确实使函数不完全是 "peak",更像是 "sharp Lorentzian minus broader Gaussian" 的混合。您可能想在 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。这是否是您真正想要做的可能取决于数据的性质。