我有一组坐标 (x, y, z(x, y))，它们描述坐标 x, y 处的强度 (z)。对于不同坐标处的一组强度，我需要拟合一个二维高斯函数来最小化均方误差。 数据位于 numpy 矩阵中，对于每个拟合会话，我将有 4、9、16 或 25 个坐标。最终我只需要得到具有最小MSE的高斯(x_0,y_0)的中心位置。 我发现的所有示例都使用 scipy.optimize.curve_fit 但它们的输入数据是整个网格而不是几个坐标。 任何帮助，将不胜感激。
有多种方法可以解决这个问题。您可以使用非线性方法（例如scipy.optimize.curve_fit），但它们会很慢并且不能保证收敛。您可以将问题线性化（快速、独特的解决方案），但分布“尾部”中的任何噪音都会导致问题。实际上，您可以在这种特殊情况下应用一些技巧来避免后一个问题。我将展示一些示例，但我现在没有时间演示所有“技巧”。
顺便说一句，一般的 2D 高斯有 6 个参数，因此您无法用 4 个点完全拟合事物。然而，听起来您可能假设 x 和 y 之间没有协方差，并且每个方向上的方差相同（即完美的“圆形”钟形曲线）。如果是这样的话，那么您只需要四个参数。如果您知道高斯的振幅，则只需要三个。不过，我将从通用解决方案开始，如果您愿意，您可以稍后对其进行简化。
目前，让我们专注于使用非线性方法（例如scipy.optimize.curve_fit).
二维高斯的一般方程是（直接来自维基百科）：
where:
is essentially 0.5 over the covariance matrix, A is the amplitude,
and (X₀, Y₀) is the center
我们把上面的等式写出来：
import numpy as np
import matplotlib.pyplot as plt
def gauss2d(x, y, amp, x0, y0, a, b, c):
inner = a * (x - x0)**2
inner += 2 * b * (x - x0)**2 * (y - y0)**2
inner += c * (y - y0)**2
return amp * np.exp(-inner)
然后让我们生成一些示例数据。首先，我们将生成一些易于拟合的数据：
np.random.seed(1977) # For consistency
x, y = np.random.random((2, 10))
x0, y0 = 0.3, 0.7
amp, a, b, c = 1, 2, 3, 4
zobs = gauss2d(x, y, amp, x0, y0, a, b, c)
fig, ax = plt.subplots()
scat = ax.scatter(x, y, c=zobs, s=200)
fig.colorbar(scat)
plt.show()
请注意，我们没有添加任何噪声，并且分布的中心位于我们拥有的数据范围内（即中心位于 0.3、0.7 处，x、y 观测值的散点在 0 和 1 之间）。目前，让我们坚持这一点，然后我们将看看当我们添加噪声并移动中心时会发生什么。
首先，让我们使用scpy.optimize.curve_fit对高斯函数进行非线性最小二乘拟合。 （顺便说一句，您可以通过使用中的一些其他函数来尝试精确的最小化算法scipy.optimize.)
The scipy.optimize函数期望的函数签名与我们上面最初编写的函数签名略有不同。我们可以编写一个包装器来“翻译”，但让我们重写一下gauss2d函数代替：
def gauss2d(xy, amp, x0, y0, a, b, c):
x, y = xy
inner = a * (x - x0)**2
inner += 2 * b * (x - x0)**2 * (y - y0)**2
inner += c * (y - y0)**2
return amp * np.exp(-inner)
我们所做的就是让函数将自变量 (x & y) 视为单个 2xN 数组。
现在我们需要初步猜测高斯曲线的参数实际上是什么。这是可选的（默认值是全一，如果我没记错的话），但是如果 1, 1 不是特别接近高斯曲线的“真实”中心，则可能会出现收敛问题。因此，我们将使用最大观测到的 z 值的 x 和 y 值作为中心的起点。我将其余参数保留为 1，但如果您知道它们可能始终存在显着差异，请将它们更改为更合理的值。
这是完整的独立示例：
import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt
def main():
x0, y0 = 0.3, 0.7
amp, a, b, c = 1, 2, 3, 4
true_params = [amp, x0, y0, a, b, c]
xy, zobs = generate_example_data(10, true_params)
x, y = xy
i = zobs.argmax()
guess = [1, x[i], y[i], 1, 1, 1]
pred_params, uncert_cov = opt.curve_fit(gauss2d, xy, zobs, p0=guess)
zpred = gauss2d(xy, *pred_params)
print 'True parameters: ', true_params
print 'Predicted params:', pred_params
print 'Residual, RMS(obs - pred):', np.sqrt(np.mean((zobs - zpred)**2))
plot(xy, zobs, pred_params)
plt.show()
def gauss2d(xy, amp, x0, y0, a, b, c):
x, y = xy
inner = a * (x - x0)**2
inner += 2 * b * (x - x0)**2 * (y - y0)**2
inner += c * (y - y0)**2
return amp * np.exp(-inner)
def generate_example_data(num, params):
np.random.seed(1977) # For consistency
xy = np.random.random((2, num))
zobs = gauss2d(xy, *params)
return xy, zobs
def plot(xy, zobs, pred_params):
x, y = xy
yi, xi = np.mgrid[:1:30j, -.2:1.2:30j]
xyi = np.vstack([xi.ravel(), yi.ravel()])
zpred = gauss2d(xyi, *pred_params)
zpred.shape = xi.shape
fig, ax = plt.subplots()
ax.scatter(x, y, c=zobs, s=200, vmin=zpred.min(), vmax=zpred.max())
im = ax.imshow(zpred, extent=[xi.min(), xi.max(), yi.max(), yi.min()],
aspect='auto')
fig.colorbar(im)
ax.invert_yaxis()
return fig
main()
在这种情况下，我们完全恢复了原始的“真实”参数。
True parameters: [1, 0.3, 0.7, 2, 3, 4]
Predicted params: [ 1. 0.3 0.7 2. 3. 4. ]
Residual, RMS(obs - pred): 1.01560615193e-16
正如我们稍后将看到的，情况并非总是如此......
让我们为我们的观察添加一些噪音。我在这里所做的就是改变generate_example_data功能：
def generate_example_data(num, params):
np.random.seed(1977) # For consistency
xy = np.random.random((2, num))
noise = np.random.normal(0, 0.3, num)
zobs = gauss2d(xy, *params) + noise
return xy, zobs
然而，结果看起来却截然不同：
就参数而言：
True parameters: [1, 0.3, 0.7, 2, 3, 4]
Predicted params: [ 1.129 0.263 0.750 1.280 32.333 10.103 ]
Residual, RMS(obs - pred): 0.152444640098
预测的中心没有太大变化，但是b and c参数改变了不少。
如果我们将函数的中心更改为稍微位于分散点之外的某个位置：
x0, y0 = -0.3, 1.1
由于存在噪音，我们最终会完全胡言乱语！ （它仍然可以正常工作，没有噪音。）
True parameters: [1, -0.3, 1.1, 2, 3, 4]
Predicted params: [ 0.546 -0.939 0.857 -0.488 44.069 -4.136]
Residual, RMS(obs - pred): 0.235664449826
这是拟合衰减为零的函数时的常见问题。 “尾巴”中的任何噪音都可能导致非常糟糕的结果。有多种策略可以解决这个问题。最简单的方法之一是通过观察到的 z 值对反演进行加权。这是一维情况的示例：（重点关注线性化问题）如何快速对多个数据集执行最小二乘拟合？ https://stackoverflow.com/questions/8780912/how-can-i-perform-a-least-squares-fitting-over-multiple-data-sets-fast/8783634#8783634如果以后有时间，我会为 2D 案例添加一个示例。