在python中用高斯拟合两个峰答案

【问题标题】：Fitting two peaks with gauss in python在python中用高斯拟合两个峰
【发布时间】：2020-07-15 22:31:57
【问题描述】：

Curve_fit 拟合不正确。我正在尝试用curve_fit 拟合实验数据。数据从.txt 文件导入到数组中：

d = np.loadtxt("data.txt")
data_x          = np.array(d[:, 0])
data_y          = np.array(d[:, 2])
data_y_err      = np.array(d[:, 3])

因为我知道必须有两个峰值，所以我的模型是两条高斯曲线的总和：

def model_dGauss(x, xc, A, y0, w, dx):
    P = A/(w*np.sqrt(2*np.pi))
    mu1 = (x - (xc - dx/3))/(2*w**2)
    mu2 = (x - (xc + 2*dx/3))/(2*w**2)    
    return y0 + P * ( np.exp(-mu1**2) + 0.5 * np.exp(-mu2**2))

使用猜测值对我的猜测值非常敏感。如果几乎完美的拟合参数可以提供结果，那么拟合数据的意义何在？还是我做错了什么？

t = np.linspace(8.4, 10, 300)
guess_dG = [32, 1, 10, 0.1, 0.2]
popt, pcov = curve_fit(model_dGauss, data_x, data_y, p0=guess_dG, sigma=data_y_err, absolute_sigma=True)
A, xc, y0, w, dx = popt

绘制数据

plt.scatter(data_x, data_y)
plt.plot(t, model_dGauss(t1,*popt))
plt.errorbar(data_x, data_y, yerr=data_y_err)

产量：

Plot result

结果只是图表底部的一条直线，而评估的参数并没有那么糟糕。怎么可能？

【问题讨论】：

sn-ps 有点不完整。尽管您应该提供数据，但您使用 t 和 t1 没有定义后者，而第一个的范围完全错误。此外，起始值似乎不太好。当放置适当的范围时，拟合效果很好。

标签： python numpy scipy curve-fitting gaussian

【解决方案1】：

一个完整的代码示例总是很受欢迎（而且，咳咳，通常在 SO 上预期）。为了消除在这里使用curve_fit 的大部分困惑，请允许我建议您使用lmfit (https://lmfit.github.io/lmfit-py) 会更轻松，尤其是它的内置模型函数和命名参数的使用。使用 lmfit，您的两个高斯加上一个恒定偏移量的代码可能如下所示：

from lmfit.models import GaussianModel, ConstantModel

# start with 1 Gaussian + Constant offset:
model = GaussianModel(prefix='p1_') + ConstantModel()

# this model will have parameters named:  
# p1_amplitude, p1_center, p1_sigma, and c.  
# here we give initial values to these parameters
params = model.make_params(p1_amplitude=10, p1_center=32, p1_sigma=0.5, c=10)

# optionally place bounds on parameters (probably not needed here):
params['p1_amplitude'].min = 0.
## params['p1_center'].vary = False # fix a parameter from varying in fit

# now do the fit (including weighting residual by 1/y_err):
result = model.fit(data_y, params, x=data_x, weights=1.0/data_y_err)

# print out param values, uncertainties, and fit statistics, or get best-fit
# parameters from `result.params`
print(result.fit_report())

# plot results
plt.errorbar(data_x, data_y, yerr=data_y_err, label='data')
plt.plot(data_x, result.best_fit, label='best fit')
plt.legend()
plt.show()

要添加第二个高斯，您可以这样做

model = GaussianModel(prefix='p1_') + GaussianModel(prefix='p2_') + ConstantModel()

# and then:  
params = model.make_params(p1_amplitude=10, p1_center=32, p1_sigma=0.5, c=10, 
                           p2_amplitude=2, p2_center=31.75, p1_sigma=0.5)

等等。

您的模型具有两个高斯共享或至少具有“链接”值 - 两个峰值的 sigma 值应该相同，并且第二个峰值的幅度是第一个峰值的一半。正如到目前为止所定义的，2-Gaussian 模型的所有参数都是独立的。但是 lmfit 有一种机制，可以通过根据其他参数给出代数表达式来对任何参数设置约束。所以，例如，你可以说

params['p2_sigma'].expr = 'p1_sigma'
params['p2_amplitude'].expr = 'p1_amplitude / 2.0'

现在，p2_amplitude 和 p2_sigma 将不会在拟合中独立变化，但会被限制为具有这些值。

【讨论】：