高斯双峰拟合

python基础

这里展示如何使用curve_fit进行高斯双峰拟合。

Pasted-image-20240905154219.png-394cf4e9c8.png

其实双峰或者多峰拟合的思路非常简单，就是把单个 gaussian 函数加起来，形成一个新的函数，再使用 curve_fit 进行拟合。比较繁琐的是，参数数量倍增。

# 定义高斯函数
def gaussian(x, a, x0, sigma, b):
    '''
    x: 输入变量
    a: 参数，缩放量
    x0：高斯峰所在位置，相当于高斯分布平均值
    sigma：高斯分布一个标准差
    b: 参数，平移量
    '''
    return a * np.exp(-(x - x0) ** 2 / (2 * sigma ** 2))+b




def double_gaussian(x, a1, mean1, sigma1, b1, a2, mean2, sigma2, b2):
    return gaussian(x, a1, mean1, sigma1, b1) + gaussian(x, a2, mean2, sigma2, b2)

然后一个很容易理解的情况是，参数越多，拟合难度越大。因此可以在 curve_fit 中设置 最大迭代次数 maxfev 这个参数为更大值（默认为1800，可以设置到20000）。

此外还可以设置拟合参数的上下限（bounds）。但如果设置了 bounds，就不能使用默认的 lm 方法进行拟合，只能用 rtf 或者 dogbox，这些在 curve_fit 的文档中都有详细说明。

然后还有一些奏效的方法，包括对数据进行预处理，减小噪声干扰；或者对高斯函数进行变形；或者更换其它更适合的算法。

适用于简单情况下的高斯双峰拟合函数代码如下：

def double_gaussian_fit(xx, yy, initial_guess, maxfev=20000):
    '''返回拟合后曲线的xy值'''
    # bounds = ((0, 0, 0, -np.inf, 0, 0, 0, -np.inf),(np.inf, 10, 10, np.inf, np.inf, 10, 10, np.inf))
    # popt, pcov = curve_fit(double_gaussian, xx, yy, p0=initial_guess, maxfev=maxfev,  method='rtf', bounds=bounds)




    popt, pcov = curve_fit(double_gaussian, xx, yy, p0=initial_guess, maxfev=maxfev)
    x3 = np.linspace(xx.min(),xx.max(), 300)  # 根据x范围均匀插值
    y3 = double_gaussian(x3, *popt)
    plt.figure(figsize=(4,3))
    plt.scatter(xx, yy, c='gray')
    plt.plot(x3, y3, c='red')
    plt.show()




    # 计算各个参数的拟合误差
    perr = np.sqrt(np.diag(pcov))




    # 计算拟合的R方值
    y2 = double_gaussian(xx, *popt)
    mean = np.mean(yy)
    ss_tot = np.sum((yy - mean) ** 2)
    ss_res = np.sum((yy - y2) ** 2)
    r_squared = 1 - (ss_res / ss_tot)
    print(f'R_squared:{r_squared}')




    loc1 = popt[1]
    loc2 = popt[5]
    sig1 = popt[2]
    sig2 = popt[6]
    FWHM1 = 2 * np.sqrt(2 * np.log(2)) * abs(sig1)
    FWHM2 = 2 * np.sqrt(2 * np.log(2)) * abs(sig2)
    # 打印拟合参数和误差
    print("参数:", popt)
    print("参数误差:", perr)
    print(f"Peak 1 location: {loc1*1000:.1f} nm\nSigma: {abs(sig1)*1000:.1f} nm\nFWHM(分辨率): {FWHM1*1000:.1f} nm")
    print(f"Peak 2 location: {loc2*1000:.1f} nm\nSigma: {abs(sig2)*1000:.1f} nm\nFWHM(分辨率): {FWHM2*1000:.1f} nm")
    return x3, y3