高斯的 2d 拟合不起作用

Question

我正在尝试拟合二维数据集。我可以毫无问题地绘制它，但是当我尝试拟合它时，我总是收到相同的错误消息，我不明白：

def twod_Gaussian(x,y,Amp,x0,y0,sigma_x,sigma_y,Offset):
    return Amp*np.exp(-((x-x0)/(2*sigma_x**2)+(y-y0)/(2*sigma_y**2)))+Offset

# get data
filename = '20180503-1455-43_confocal_xy_data.dat'
data = np.genfromtxt(filename, comments='#')
values = open(filename)
val = values.readlines(500)
res = float(val[11][-4:-1])
values.close()
x = [float(i) for i in data[:, 0]]
y = [float(i) for i in data[:, 1]]
count_rate = [float(i) for i in data[:, 3]]



#reshape data
xmax = max(x)
xmin = min(x)
ymax = max(y)
ymin = min(y)
cmax = max(count_rate)
cmin = min(count_rate)
np.percentile(count_rate, 99)
resX = (xmax-xmin)/res
rexY = (ymax-ymin)/res
res2 = len(x)/res
xx = np.array(x).reshape((int(res2), int(res)))
yy = np.array(y).reshape((int(res2), int(res)))



#plot data
count_matrix = np.array(count_rate).reshape((int(res2), int(res)))
fig = plt.figure()
ax = fig.gca(projection='3d')
CS = ax.plot_surface((xx)/1e-6, (yy)/1e-6, count_matrix, cmap='plasma')
plt.ticklabel_format(style='plain', axis='both', scilimits=(0, 0))
ax.set_xlabel(r'X position ($\mathrm{\mu}$m)')
ax.set_ylabel(r'Y position ($\mathrm{\mu}$m)')
fig.colorbar(CS, ax=ax, extend='max', format='%.0e')
plt.show()

#initial guess and trying to fit it
p0=(98e-06,81.5e-06,50000,0.5e-06,0.5e-06,20000)
popt, pcov = curve_fit(twod_Gaussian, x, y, count_rate, p0)

我总是收到错误消息

ValueError: sigma 的形状不正确。

如果我将最后一行更改为

popt, pcov = curve_fit(twod_Gaussian, (x, y), count_rate, p0)

我收到错误消息：

twod_Gaussian() missing 1 required positional argument: 'Offset'

我不知道为什么形状应该是错误的？

Answer 1

首先，我认为问题在于您的模型函数定义有一个额外的参数，因为您对二维使用两个单独的参数（x 和 y），而不是独立变量的一个参数curve_fit 需要。将您的模型函数更改为

def twod_Gaussian(xy, Amp,x0,y0,sigma_x,sigma_y,Offset):
    x, y = xy
    return Amp*np.exp(-((x-x0)/(2*sigma_x**2)+(y-y0)/(2*sigma_y**2)))+Offset

然后调用curve_fit

xy = x, y
popt, pcov = curve_fit(twod_Gaussian, xy, count_rate, p0)

应该做你想做的。

其次，您也可以考虑使用lmfit (https://lmfit.github.io/lmfit-py)。这个库有一个稍微高级一点的曲线拟合接口，仍然围绕着scipy.optimize求解器。例如，它将参数视为一等的命名对象，并轻松支持在模型函数之外对参数设置边界和约束。对于您的情况，它还支持为拟合函数定义多个自变量，而不是强迫您假设只有模型函数的第一个参数是自变量，而所有其他参数都是拟合中的变量。使用 lmfit，您的问题看起来大致如下（我没有运行它，因为我没有您的数据）：

from lmfit import Model
def twod_Gaussian(x, y, Amp,x0,y0,sigma_x,sigma_y,Offset):
    return Amp*np.exp(-((x-x0)/(2*sigma_x**2)+(y-y0)/(2*sigma_y**2)))+Offset

# turn your model function into a Model, specifying independent variables
g2model = Model(twod_Gaussian, independent_vars=['x', 'y'])

# create Parameters object -- ordered dict of *named parameters, 
# giving initial values here:
params = g2model.make_params(Amp=98e-6, x0=81e-6, y0=50000., sigma_x=5e-7, sigma_y=5e-7, Offset=20000)

# Note that using names is *way* better than using an ordered list. 
# for example, you may have mixed up the order and really meant:
params = g2model.make_params(x0=98e-6, y0=81e-6, Amp=50000., sigma_x=5e-7, sigma_y=5e-7, Offset=20000)

# you can place bounds on parameters, perhaps as
params['sigma_x'].min = 0
params['sigma_y'].min = 0

# do fit, passing independent vars explicitly
result = g2model.fit(count_rate, params, x=x, y=y)

# print out full report of statistics, best-fit values and stderrs:
print(result.fit_report())

# array of best-fit data == result.best_fit

文档中描述了更多功能。

最后，要使用您的函数制作二维高斯，y 应该是按行排列的数组（x 是按列排列的）。否则你的twod_Gaussian() 将返回一个与x 和y 长度相同的一维数组。所以，我想你想添加

y = np.vstack(y)

就在您从 data 数组创建 y 之后。