Python双高斯拟合异常：第二个高斯曲线形态不符数据

我是Python新手，尝试用数据拟合双高斯函数，能画出拟合图像，但第二个高斯曲线形态远大于数据实际情况，不符合预期。以下是实现代码及拟合结果图，求解决方法：

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import scipy as scipy
from scipy import optimize
from matplotlib.ticker import AutoMinorLocator
from matplotlib import gridspec
import matplotlib.ticker as ticker
%matplotlib inline

data = np.loadtxt('csv/test_run09.csv', encoding="utf-8", delimiter=',',skiprows=1)
x = data[:,1]
y1 = data[:,2]
y2 = data[:,3]
y3 = data[:,4]
y4 = data[:,5]
y5 = data[:,6]
y6 = data[:,7]
y7 = data[:,8]
y8 = data[:,9]
y9 = data[:,10]
y10 = data[:,11]
y11 = data[:,12]
y12 = data[:,13]
y13 = data[:,14]
y14 = data[:,15]

amp1 = 2
sigma1 = 0.1
x_array = x[(x>33)&(x<34)]
y_array = y14[(x>33)&(x<34)]

amp2 = np.max(y_array)
sigma2 = np.std(x_array)

def _1gaussian1(x_array, amp1, sigma1):
    return amp1*(1/(sigma1*(np.sqrt(2*np.pi))))*(np.exp((-1.0/2.0)*(((x_array-33.49290958)/sigma1)**2))) + 0.2
def _1gaussian2(x_array, amp2, sigma2):
    return amp2*(1/(sigma2*(np.sqrt(2*np.pi))))*(np.exp((-1.0/2.0)*(((x_array-33.6312849)/sigma2)**2))) + 0.2

popt_gauss1, pcov_gauss1 = scipy.optimize.curve_fit(_1gaussian1, x_array, y_array, p0=[amp1, sigma1])
popt_gauss2, pcov_gauss2 = scipy.optimize.curve_fit(_1gaussian2, x_array, y_array, p0=[np.max(y_array), np.std(x_array)])

def _2gaussian(x_array, amp1, sigma1, amp2, sigma2):
    return amp1*(1/(sigma1*(np.sqrt(2*np.pi))))*(np.exp((-1.0/2.0)*(((x_array-33.49290958)/sigma1)**2))) + amp2*(1/(sigma1*(np.sqrt(2*np.pi))))*(np.exp((-1.0/2.0)*(((x_array-33.6312849)/sigma2)**2))) + 0.3

popt_2gauss, pcov_2gauss = scipy.optimize.curve_fit(_2gaussian, x_array, y_array, p0=[amp1, sigma1, np.max(y_array), np.std(x_array)])
perr_2gauss = np.sqrt(np.diag(pcov_2gauss))
print(popt_2gauss)

pars_1 = popt_2gauss[0:2]
pars_2 = popt_2gauss[2:4]
gauss_peak_1 = _1gaussian1(x_array, *pars_1)
gauss_peak_2 = _1gaussian2(x_array, *pars_2)

fig = plt.figure(figsize=(7,5))
gs = gridspec.GridSpec(1,1)
ax1 = fig.add_subplot(gs[0])
plt.grid()

ax1.plot(x_array, y_array, "ro")
ax1.plot(x_array, _2gaussian(x_array, *popt_2gauss), 'k--')

# peak 1
ax1.plot(x_array, gauss_peak_1, "g")
ax1.fill_between(x_array, gauss_peak_1.min(), gauss_peak_1, facecolor="green", alpha=0.5)

# peak 2
ax1.plot(x_array, gauss_peak_2, "y")
ax1.fill_between(x_array, gauss_peak_2.min(), gauss_peak_2, facecolor="yellow", alpha=0.5)  

# prints the fitting parameters with their errors
print("-------------Peak 1-------------")
print("amplitude = %0.2f (+/-) %0.2f" % (pars_1[0], perr_2gauss[0]))
print("sigma = %0.2f (+/-) %0.2f" % (pars_1[1], perr_2gauss[1]))
print("area = %0.2f" % np.trapz(gauss_peak_1))
print("-------------Peak 2-------------")
print("amplitude = %0.2f (+/-) %0.2f" % (pars_2[0], perr_2gauss[2]))
print("sigma = %0.2f (+/-) %0.2f" % (pars_2[1], perr_2gauss[3]))
print("area = %0.2f" % np.trapz(gauss_peak_2))

问题解决步骤 #

1. 修复双高斯函数的公式错误
   你的_2gaussian函数里，第二个高斯项的分母错误使用了sigma1，应该替换为sigma2：

   def _2gaussian(x_array, amp1, sigma1, amp2, sigma2):
       return amp1*(1/(sigma1*(np.sqrt(2*np.pi))))*np.exp(-0.5*((x_array-33.49290958)/sigma1)**2) + \
              amp2*(1/(sigma2*(np.sqrt(2*np.pi))))*np.exp(-0.5*((x_array-33.6312849)/sigma2)**2) + 0.3
   

2. 修正单高斯拟合的调用错误
   代码里scipy.optimize.curve_fit(_2gaussian2, ...)写错了函数名，应该是_1gaussian2，这个错误会导致无法正确得到第二个单高斯的拟合参数，进而影响双高斯的初始值。

3. 优化初始参数设置
   不要直接用np.max(y_array)和np.std(x_array)作为第二个高斯的初始参数，应该用单高斯拟合得到的结果作为双高斯的初始值：

   p0 = [popt_gauss1[0], popt_gauss1[1], popt_gauss2[0], popt_gauss2[1]]
   

4. 统一基线设置
   单高斯函数里加了0.2的基线，双高斯里加了0.3，这会导致拟合时基线不一致。建议把基线也作为拟合参数，或者统一为同一个值。比如修改双高斯函数，把基线设为可拟合参数：

   def _2gaussian(x_array, amp1, sigma1, amp2, sigma2, baseline):
       return amp1*(1/(sigma1*(np.sqrt(2*np.pi))))*np.exp(-0.5*((x_array-33.49290958)/sigma1)**2) + \
              amp2*(1/(sigma2*(np.sqrt(2*np.pi))))*np.exp(-0.5*((x_array-33.6312849)/sigma2)**2) + baseline
   

   然后初始参数加上基线的初始值（比如0.2）：

   p0 = [popt_gauss1[0], popt_gauss1[1], popt_gauss2[0], popt_gauss2[1], 0.2]
   

5. 约束参数范围
   使用curve_fit的bounds参数限制参数的合理范围，比如sigma必须大于0，振幅不能为负：

   bounds = (0, [np.inf, 0.5, np.inf, 0.5, 1])  # 对应amp1, sigma1, amp2, sigma2, baseline的上下限
   popt_2gauss, pcov_2gauss = scipy.optimize.curve_fit(_2gaussian, x_array, y_array, p0=p0, bounds=bounds)
   

内容的提问来源于stack exchange，提问作者EL-san