python - 2D高斯拟合在Python中某些坐标处的强度

Question

我有一组坐标(x，y，z(x，y))，它们描述了在坐标x，y处的强度(z)。对于在不同坐标下的一定数量的这些强度，我需要拟合一个使均方误差最小的2D高斯函数。
数据以numpy矩阵表示，对于每个拟合 session ，我将具有4、9、16或25个坐标。最终，我只需要获得MSE最小的高斯(x_0，y_0)的中心位置。
我发现的所有示例都使用scipy.optimize.curve_fit，但是它们拥有的输入数据是在整个网格上而不是几个坐标上。
任何帮助，将不胜感激。

Answer 1

介绍

有多种方法可以解决此问题。您可以使用非线性方法(例如scipy.optimize.curve_fit)，但是它们会很慢并且不能保证收敛。您可以线性化问题(快速，独特的解决方案)，但是分布“尾部”中的任何噪声都会引起问题。实际上，您可以将一些技巧应用于此特定情况，以避免出现后一个问题。我将展示一些示例，但是我现在没有时间展示所有的“技巧”。

顺便提一句，一般的2D波斯语有6个参数，因此您将无法完全用4个点拟合事物。但是，听起来您可能假设x和y之间没有协方差，并且每个方向的方差都相同(即完美的“圆”钟形曲线)。如果是这样，那么您只需要四个参数。如果您知道高斯的幅度，则只需三个即可。但是，我将从通用解决方案开始，如果需要，您可以稍后对其进行简化。

目前，让我们着重于使用非线性方法(例如scipy.optimize.curve_fit)解决此问题。

二维高斯的一般方程为(直接来自维基百科):

在哪里:

在协方差矩阵上基本上是0.5，A是振幅，
并且(X₀，Y₀)是中心

生成简化的样本数据

让我们将上面的方程写出来:

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中的其他一些功能来尝试使用精确的最小化算法。)
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”参数。

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和 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值对反演进行加权。这是一维情况的示例:(专注于线性化问题) How can I perform a least-squares fitting over multiple data sets fast?如果有时间，我将为2D情况添加一个示例。