0

我正在研究一些需要能够执行二维高斯拟合的代码。我的代码主要基于以下问题：Fitting a 2D Gaussian function using scipy.optimize.curve_fit - ValueError and minpack.error。现在的问题是我对需要使用的不同参数并没有真正的初步猜测。

我试过这个：

def twoD_Gaussian(x_data_tuple, amplitude, xo, yo, sigma_x, sigma_y, theta, offset):
    (x,y) = x_data_tuple
    xo = float(xo)
    yo = float(yo)    
    a = (np.cos(theta)**2)/(2*sigma_x**2) + (np.sin(theta)**2)/(2*sigma_y**2)
    b = -(np.sin(2*theta))/(4*sigma_x**2) + (np.sin(2*theta))/(4*sigma_y**2)
    c = (np.sin(theta)**2)/(2*sigma_x**2) + (np.cos(theta)**2)/(2*sigma_y**2)
    g = offset + amplitude*np.exp( - (a*((x-xo)**2) + 2*b*(x-xo)*(y-yo) 
                            + c*((y-yo)**2)))
    return g.ravel()

data.reshape(201,201) 只是我从上述问题中得到的。

mean_gauss_x = sum(x * data.reshape(201,201)) / sum(data.reshape(201,201))
sigma_gauss_x = np.sqrt(sum(data.reshape(201,201) * (x - mean_gauss_x)**2) / sum(data.reshape(201,201)))

mean_gauss_y = sum(y * data.reshape(201,201)) / sum(data.reshape(201,201))
sigma_gauss_y = np.sqrt(sum(data.reshape(201,201) * (y - mean_gauss_y)**2) / sum(data.reshape(201,201)))


initial_guess = (np.max(data), mean_gauss_x, mean_gauss_y, sigma_gauss_x, sigma_gauss_y,0,10)


popt, pcov = curve_fit(twoD_Gaussian, (x, y), data, p0=initial_guess)

data_fitted = twoD_Gaussian((x, y), *popt)

如果我尝试这样做，我会收到以下错误消息：ValueError: setting an array element with a sequence。

关于开始参数的推理是否正确？为什么我会收到这个错误？

4

1 回答 1

1

如果我使用链接问题中的可运行代码并替换您的定义initial_guess：

mean_gauss_x = sum(x * data.reshape(201,201)) / sum(data.reshape(201,201))
sigma_gauss_x = np.sqrt(sum(data.reshape(201,201) * (x - mean_gauss_x)**2) / sum(data.reshape(201,201)))

mean_gauss_y = sum(y * data.reshape(201,201)) / sum(data.reshape(201,201))
sigma_gauss_y = np.sqrt(sum(data.reshape(201,201) * (y - mean_gauss_y)**2) / sum(data.reshape(201,201)))

initial_guess = (np.max(data), mean_gauss_x, mean_gauss_y, sigma_gauss_x, sigma_gauss_y,0,10)

然后

print(inital_guess)

产量

(13.0, array([...]), array([...]), array([...]), array([...]), 0, 10)

请注意，其中的一些值initial_guess是数组。该optimize.curve_fit函数期望initial_guess是一个标量元组。这就是问题的根源。

错误信息

ValueError: setting an array element with a sequence

当预期标量值时提供类似数组时，通常会出现这种情况。这暗示问题的根源可能与维数错误的数组有关。例如，如果您将一维数组传递给需要标量的函数，则可能会出现这种情况。

让我们看一下取自链接问题的这段代码：

x = np.linspace(0, 200, 201)
y = np.linspace(0, 200, 201)
X, Y = np.meshgrid(x, y)

x和y是一维数组，而X和Y是二维数组。（我已将所有二维数组大写以帮助将它们与一维数组区分开来）。

现在请注意，Pythonsum和 NumPy 的sum方法在应用于 2D 数组时表现不同：

In [146]: sum(X)
Out[146]: 
array([    0.,   201.,   402.,   603.,   804.,  1005.,  1206.,  1407.,
        1608.,  1809.,  2010.,  2211.,  2412.,  2613.,  2814.,  3015.,
        ...
       38592., 38793., 38994., 39195., 39396., 39597., 39798., 39999.,
       40200.])

In [147]: X.sum()
Out[147]: 4040100.0

Pythonsum函数相当于

total = 0
for item in X:
    total += item

由于X是 2D 数组，因此循环for item in X遍历X 的行。因此，每个item都是代表一行的一维数组X。因此，total最终成为一维数组。

相反，X.sum()将所有元素相加X并返回一个标量。

因为initial_guess应该是一个标量元组，所以在你使用的任何地方都sum应该使用 NumPysum方法。例如，替换

mean_gauss_x = sum(x * data) / sum(data)

和

mean_gauss_x = (X * DATA).sum() / (DATA.sum())

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

# define model function and pass independant variables x and y as a list
def twoD_Gaussian(data, amplitude, xo, yo, sigma_x, sigma_y, theta, offset):
    X, Y = data
    xo = float(xo)
    yo = float(yo)
    a = (np.cos(theta) ** 2) / (2 * sigma_x ** 2) + (np.sin(theta) ** 2) / (
        2 * sigma_y ** 2
    )
    b = -(np.sin(2 * theta)) / (4 * sigma_x ** 2) + (np.sin(2 * theta)) / (
        4 * sigma_y ** 2
    )
    c = (np.sin(theta) ** 2) / (2 * sigma_x ** 2) + (np.cos(theta) ** 2) / (
        2 * sigma_y ** 2
    )
    g = offset + amplitude * np.exp(
        -(a * ((X - xo) ** 2) + 2 * b * (X - xo) * (Y - yo) + c * ((Y - yo) ** 2))
    )
    return g.ravel()


# Create x and y indices
x = np.linspace(0, 200, 201)
y = np.linspace(0, 200, 201)
X, Y = np.meshgrid(x, y)

# create data
data = twoD_Gaussian((X, Y), 3, 100, 100, 20, 40, 0, 10)
data_noisy = data + 0.2 * np.random.normal(size=data.shape)
DATA = data.reshape(201, 201)


# add some noise to the data and try to fit the data generated beforehand
mean_gauss_x = (X * DATA).sum() / (DATA.sum())
sigma_gauss_x = np.sqrt((DATA * (X - mean_gauss_x) ** 2).sum() / (DATA.sum()))

mean_gauss_y = (Y * DATA).sum() / (DATA.sum())
sigma_gauss_y = np.sqrt((DATA * (Y - mean_gauss_y) ** 2).sum() / (DATA.sum()))


initial_guess = (
    np.max(data),
    mean_gauss_x,
    mean_gauss_y,
    sigma_gauss_x,
    sigma_gauss_y,
    0,
    10,
)
print(initial_guess)
# (13.0, 100.00000000000001, 100.00000000000001, 57.106515650488404, 57.43620227324201, 0, 10)
# initial_guess = (3,100,100,20,40,0,10)

popt, pcov = optimize.curve_fit(twoD_Gaussian, (X, Y), data_noisy, p0=initial_guess)

data_fitted = twoD_Gaussian((X, Y), *popt)

fig, ax = plt.subplots(1, 1)
ax.imshow(
    data_noisy.reshape(201, 201),
    cmap=plt.cm.jet,
    origin="bottom",
    extent=(X.min(), X.max(), Y.min(), Y.max()),
)
ax.contour(X, Y, data_fitted.reshape(201, 201), 8, colors="w")
plt.show()
于 2019-08-13T15:19:15.383 回答