使用一维方法进行 2D scipy.optimize.curve_fit

Question

我的目标是使用 scipy curve_fit 函数实现二维曲线拟合（我有点懒，只是想将典型的 1D 曲线拟合应用于 2D 曲面）。

我从两个 lambda 函数开始。第一个描述了带有变量 X、Y 和参数 P 的二维表面； lambda1(X,Y,P)。第二个 lambda 函数“包装”第一个 lambda，这样我就可以在一维数组上使用它，lambda2(index, P) = lambda1(x[index], y[index], P)。然后，我编写 3 个二维矩阵：XX、YY、ZZ 矩阵，并将其重塑为三个一维数组：x_line、y_line、z_line。

lambda 函数似乎表现良好，因为当我调用 lambda2 并将其输出重塑为正确的形式时，我可以创建一个二维表面。然而，当我调用 curve_fit 来查找最佳参数时，我遇到了错误提示： ”
sag_eq_wrapped = lambda 索引, x0, y0, z0, r, k: ((1/r)((x_line[index]-x0)**2+(y_line[index]-y0)**2)**2 ) / (1+np.sqrt(1-(1+k)(1/r)*2((x_line[索引]-x0)**2+(y_line[索引]-y0)**2 ))) + z0 ~~~~~~^^^^^^^ IndexError：用作索引的数组必须是整数（或布尔）类型 ”

有办法解决这个错误吗？我之前在 curve_fit 中使用过 lambda 函数，但这是我第一次尝试这种包装方法并在 lambda 函数中调用其他数组。也许它只是与 curve_fit 函数表现不佳......

感谢您的帮助。

下面是示例代码：

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

###Define Lanbda function for Surface
sag_eq = lambda x, y, x0, y0, z0, r, k: ((1/r)*((x-x0)**2+(y-y0)**2)**2) / (1+np.sqrt(1-(1+k)*(1/r)**2*((x-x0)**2+(y-y0)**2))) + z0
sag_eq_wrapped = lambda index, x0, y0, z0, r, k: ((1/r)*((x_line[index]-x0)**2+(y_line[index]-y0)**2)**2) / (1+np.sqrt(1-(1+k)*(1/r)**2*((x_line[index]-x0)**2+(y_line[index]-y0)**2))) + z0

###Make Image and Domain
x_dom = np.linspace(0,0.32, 2**8)
y_dom = np.linspace(0,0.24, 2**7)
param = [x_dom[-1]/2, y_dom[-1]/2, 0, -0.2, -1]
zz = []
for y in y_dom:
    base = sag_eq(x_dom, y, param[0], param[1], param[2], param[3], param[4])
    noise = np.random.rand(len(base))
    zz.append(base+noise*1e-3)

###Make the two-dimnesional Data
zz=np.array(zz)
xx, yy = np.meshgrid(x_dom, y_dom)

###Make wrapped 1Dimnesional arrays
x_line = xx.reshape(len(x_dom)*len(y_dom))
y_line = yy.reshape(len(x_dom)*len(y_dom))
z_line = zz.reshape(len(x_dom)*len(y_dom))
index_line = np.arange(0, len(x_dom)*len(y_dom), dtype=int)

###Make surface using wrapped Lambda function
surface_1 = sag_eq_wrapped(index_line, param[0], param[1], param[2], param[3], param[4])
surface_1 = surface_1.reshape(len(y_dom), len(x_dom))

fig, [ax_data, ax_surf1] = plt.subplots(1,2)
ax_data.contourf(xx, yy, zz); ax_data.set_title("Data Surf")
ax_surf1.contourf(xx, yy, surface_1); ax_surf1.set_title("lamda made Surf")
plt.show()

fig, [ax_x, ax_y] = plt.subplots(1, 2)
ax_x.scatter(x_line, z_line); ax_x.set_title("Wrapped X Data")
ax_y.scatter(y_line, z_line); ax_y.set_title("Wrapped Y Data")
plt.show()


popt, pcov = curve_fit(sag_eq_wrapped, index_line, z_line, p0 = param)

我期望 Curve_fit 能够正常工作。不幸的是我不完全理解错误提示，因为index_line中的元素是int类型

Answer 1

方法

curve_fit

可以拟合多个自变量函数。您所需要的只是实现目标函数签名。

给定的模型具有正确的签名：

def model(x, x0, y0, z0, r, k):
    return ((1 / r) * ((x[0] - x0) ** 2 + (x[1] - y0) **2) **2) / (1 + np.sqrt(1 - (1 + k) * (1/r) ** 2 * ((x[0] - x0) ** 2 + (x[1] - y0) ** 2))) + z0

我们可以用它来适应您的数据集：

popt, pcov = optimize.curve_fit(model, [x_line, y_line], z_line, p0=[0.2, 0.2, 0., -0.1, -1])
#(array([ 1.59933408e-01,  1.20135077e-01,  4.99474227e-04, -1.98327589e-01,
#         -1.05312362e+00]),
# array([[ 6.23363842e-09, -9.67064719e-14, -9.17427476e-14,
#         -2.12181236e-10,  7.81917668e-09],
#        [-9.67064719e-14,  1.40063897e-08,  2.09623759e-13,
#          2.83771229e-10, -1.24165280e-08],
#        [-9.17427476e-14,  2.09623759e-13,  7.59814526e-12,
#          4.16322029e-09, -9.09518272e-08],
#        [-2.12181236e-10,  2.83771229e-10,  4.16322029e-09,
#          4.35907326e-06, -1.07458297e-04],
#        [ 7.81917668e-09, -1.24165280e-08, -9.09518272e-08,
#         -1.07458297e-04,  2.77570575e-03]]))

并回归函数：

zhat = model([x_line, y_line], *popt)

从图形上看，它导致：

fig, axe = plt.subplots(subplot_kw={"projection": "3d"})
axe.scatter(x_line, y_line, z_line, marker=".", alpha=0.15)
axe.plot_surface(xx, yy, zhat.reshape(xx.shape), cmap="jet")