我有各种光谱图，我正在尝试将函数拟合到其中。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。
#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-度。这是否是你真正想做的，可能取决于你的数据的性质。