我很难理解为什么我的高斯拟合一组数据(ydata
)如果我移动对应于该数据的x值的间隔(xdata1
到xdata2
)不能很好地工作。高斯写成:
其中A只是一个幅度因子。改变数据的某些值,很容易使它适用于这两种情况,但是人们也可以很容易地找到它对xdata1
不能很好地工作的情况,以及不估计参数的协方差。我在Windows 7上使用Spyder中的scipy.optimize.curve_fit
和Python 3.7.1。
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
xdata1 = np.linspace(-9,4,20, endpoint=True) # works fine
xdata2 = xdata1+2
ydata = np.array([8,9,15,12,14,20,24,40,54,94,160,290,400,420,300,130,40,10,8,4])
def gaussian(x, amp, mean, sigma):
return amp*np.exp(-(((x-mean)**2)/(2*sigma**2)))/(sigma*np.sqrt(2*np.pi))
popt1, pcov1 = curve_fit(gaussian, xdata1, ydata)
popt2, pcov2 = curve_fit(gaussian, xdata2, ydata)
fig, ([ax1, ax2]) = plt.subplots(nrows=1, ncols=2,figsize=(9, 4))
ax1.plot(xdata1, ydata, 'b+:', label='xdata1')
ax1.plot(xdata1, gaussian(xdata1, *popt1), 'r-', label='fit')
ax1.legend()
ax2.plot(xdata2, ydata, 'b+:', label='xdata2')
ax2.plot(xdata2, gaussian(xdata2, *popt2), 'r-', label='fit')
ax2.legend()
问题是你在拟合高斯时的第二次尝试是在搜索参数空间时陷入局部最小值:curve_fit是least_squares的包装器,它使用梯度下降来最小化成本函数,这很容易得到stuck in local minima。
您应该尝试提供合理的起始参数(通过使用p0
的curve_fit参数)来避免这种情况:
#... your code
y_max = np.max(y_data)
max_pos = ydata[ydata==y_max][0]
initial_guess = [y_max, max_pos, 1] # amplitude, mean, std
popt2, pcov2 = curve_fit(gaussian, xdata2, ydata, p0=initial_guess)
你可以看到哪一个合理的适合:
您应该编写一个可以提供合理的起始参数估计的函数。在这里,我刚刚找到了最大y值,并用它来确定初始参数。我发现这适用于拟合正态分布但你可以考虑其他方法。
编辑:
您也可以通过缩放幅度来解决问题:幅度太大,参数空间失真,梯度下降只是跟随幅度的最大变化方向,并有效地忽略了西格玛。请考虑参数空间中的以下图表(颜色是给定参数拟合的平方残差之和,白色十字表示最佳解决方案):
确保记下x和y轴的不同比例。
需要在y(幅度)上制作大量“单位”大小的步骤,以便从点x,y =(0,0)得到最小值,其中您只需要少于一个“单位”大小的步骤达到x(sigma)的最小值。该算法简单地采用幅度步长,因为这是最陡的梯度。当它达到最小化成本函数的幅度时,它只是停止算法,因为它似乎已经收敛并且在sigma参数中几乎没有或没有变化。
解决此问题的一种方法是缩放你的ydata以使参数空间失真:将你的ydata
除以100,你会发现你的合体无需提供任何启动参数!