使用scipy curve_fit做高斯拟合时21点部分数据集拟合错误问题

问题背景 #

• 使用scipy的curve_fit函数对数据集进行高斯函数拟合
• 测试数据集分为两类：第一类包含19个数据点，第二类包含21个数据点，两类数据集均覆盖0.5-0.7、1.0-1.2、1.5-1.7三个区间的子数据集
• 拟合异常现象：19点数据集的三个区间子数据集全部拟合成功，21点数据集仅1.5-1.7区间拟合结果正确，其余区间拟合偏差极大

原始拟合代码 #

import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

#function declaration
def gauss(x, amp, mu, sigma):
    y = amp*np.exp(-(x-mu)**2/(2*sigma**2))
    return y

#fitting
popt, pcov = curve_fit(f = gauss, xdata = x, ydata = y)
#print(popt)
amp = popt[0]
mu = popt[1]
sigma = popt[2]
print(amp,mu,sigma)

#krypton value
krypton_y = amp/((math.exp(1))**2)
#print(krypton_y)
krypton_x1 = mu + math.sqrt((-2*(sigma**2))*math.log(krypton_y/amp))
krypton_x2 = mu - math.sqrt((-2*(sigma**2))*math.log(krypton_y/amp))
print(krypton_x1-krypton_x2)

#print(gauss([krypton_x1, krypton_x2], popt[0], popt[1], popt[2]))


#horizontal line
horizontal_x = np.arange(min(x)-0.01, max(x)+0.02, 0.01)
horizontal_y = np.repeat(0, len(horizontal_x))


#build fit set
x_test = np.arange(min(x), max(x), 0.0000001)
y_test = gauss(x_test, popt[0], popt[1], popt[2])
y_krypton = []
for i in horizontal_x:
    y_krypton.append(krypton_y)


#Vertical lines
vertical_y = np.arange(-20, amp+20, 0.01)
l = len(vertical_y)
vertical_mean = np.repeat(mu, l)


#fit data
fig = plt.figure()
fig = plt.scatter(x,y, label ='original data', color = 'red', marker = 'x')
fig = plt.plot(x_test, y_test, label = 'Gaussian fit curve')
fig = plt.plot(horizontal_x, y_krypton, color = '#830000', linewidth = 1) 
fig = plt.plot(vertical_mean, vertical_y, color = '#0011ed')
fig = plt.xlabel('Distance in mm')
fig = plt.ylabel('Current in nA')
fig = plt.title('Intensity Profile for '+gas+' laser | Z = '+str(z)+'cm')
fig = plt.scatter(mu, amp, s = 25, color = '#0011ed')
fig = plt.scatter(krypton_x1, krypton_y, s = 25,  color = '#830000')
fig = plt.scatter(krypton_x2, krypton_y, s = 25,  color = '#830000')
plt.annotate('('+"{:.4f}".format(mu)+','+"{:.4f}".format(amp)+')', (mu, amp), xytext = (mu+0.002,amp+0.5))
plt.annotate('('+"{:.4f}".format(krypton_x1)+','+"{:.4f}".format(krypton_y)+')', (krypton_x1, krypton_y), xytext = (krypton_x1+0.002,krypton_y+0.5))
plt.annotate('('+"{:.4f}".format(krypton_x2)+','+"{:.4f}".format(krypton_y)+')', (krypton_x2, krypton_y), xytext = (krypton_x2+0.002,krypton_y+0.5))
plt.legend()
plt.margins(0)
plt.show()

拟合结果说明 #

• 1.5-1.7区间拟合结果正常：拟合曲线完全贴合原始数据点，峰值、半高宽计算准确
• 1.0-1.2区间拟合结果异常：拟合曲线偏离原始数据点，参数计算误差极大

问题原因 #

1. curve_fit默认使用全1作为初始拟合参数，未结合数据集本身的分布特征设置初值，非线性最小二乘求解过程中容易陷入局部最优解。21点数据集0.5-0.7、1.0-1.2区间的峰值、中心位置等特征和默认初值偏差较大，因此拟合失败；1.5-1.7区间特征与默认初值匹配度高，因此拟合正常。19点数据集所有区间特征均与默认初值偏差小，因此全部拟合成功。
2. 现有高斯函数未设置基线偏移参数，若数据存在非零基线，也会放大拟合误差。

解决方案 #

1. 设置合理的初始拟合参数 #

根据数据集的统计特征计算初值，传入curve_fit的p0参数：

• 幅度amp初值取y数据的最大值：p0_amp = np.max(y)
• 中心mu初值取y最大值对应的x值：p0_mu = x[np.argmax(y)]
• 标准差sigma初值取x区间跨度的1/4：p0_sigma = (np.max(x)-np.min(x))/4
  调整后的拟合调用代码：

popt, pcov = curve_fit(f = gauss, xdata = x, ydata = y, p0=[p0_amp, p0_mu, p0_sigma])

2. 优化高斯函数，增加基线参数 #

若数据存在非零基线，可将高斯函数修改为带偏移的形式：

def gauss(x, amp, mu, sigma, offset):
    y = amp*np.exp(-(x-mu)**2/(2*sigma**2)) + offset
    return y

对应增加offset初值，取y数据的最小值即可：p0_offset = np.min(y)，调用curve_fit时将其加入p0列表。

3. 增加参数边界限制（可选） #

为避免拟合出不合理的参数值，可设置参数上下界：

# 边界格式：([参数最小值列表], [参数最大值列表])
bounds = ([0, np.min(x), 0, -np.inf], [np.inf, np.max(x), np.inf, np.inf])
popt, pcov = curve_fit(f = gauss, xdata = x, ydata = y, p0=[p0_amp, p0_mu, p0_sigma, p0_offset], bounds=bounds)

以上调整完成后，即可解决不同区间的拟合偏差问题。

内容的提问来源于stack exchange，提问作者Tathagata Dey