用 2D 高斯拟合 2D 直方图

Question

我的问题看起来很简单，但我已经坚持了一段时间了。 我只是想用 2D 高斯总和来拟合 2D 直方图。

由于似乎没有任何效果，我回到基础知识并尝试将 2D 高斯分布与 2D 高斯拟合。 我注意到，如果我拟合人工生成的噪声数据，我不会得到与采样高斯分布并将其合并到二维直方图中时相同的结果。即使初始猜测参数非常好（比我在实际应用问题中得到的要好得多），后者也会导致质量非常差。

这是示例代码：

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

def twoD_Gaussian(xy, amplitude, xo, yo, sigma_x, sigma_y, rho):
    x, y = xy
    xo = float(xo)
    yo = float(yo)
    a = 1 / sigma_x**2
    b = rho / sigma_x / sigma_y
    c = 1 / sigma_y**2
    g = amplitude * np.exp( - (1/(2*(1-rho**2)))*(a*((x-xo)**2) - 2*b*(x-xo)*(y-yo) + c*((y-yo)**2)))
    return g.ravel()

# Produce a number of points in x-y from 1 distribution. 
mean = [0,0]
cov = [[3,1.5],[1.5,1]] 
N = int(1e6)
x,y = np.random.multivariate_normal(mean,cov,N).T

nbr_bins = 100 #int(2 * N ** (1/3))
print("numbins", nbr_bins)

# Produce 2D histogram
H,xedges,yedges = np.histogram2d(x,y,bins=nbr_bins)
bin_centers_x = (xedges[:-1]+xedges[1:])/2.0
bin_centers_y = (yedges[:-1]+yedges[1:])/2.0
X,Y = np.meshgrid(bin_centers_x, bin_centers_y)

data = twoD_Gaussian((X, Y), H.max(), mean[0], mean[1], np.sqrt(cov[0][0]), np.sqrt(cov[1][1]), cov[0][1]/(np.sqrt(cov[0][0])*np.sqrt(cov[1][1])))
data_noisy = data + 0.2*np.random.normal(size=data.shape)

# Initial Guess
p0 = (H.max(),mean[0],mean[1],1.7,1,0.9)

# Curve Fit parameters with histo
coeff, var_matrix = curve_fit(twoD_Gaussian,(X,Y),H.ravel(),p0=p0)

# Curve fit on noisy data
popt, pcov = curve_fit(twoD_Gaussian, (X, Y), data_noisy, p0=p0)

print('hist fit', coeff)
print('noisy data fit', popt)

data_fitted_hist = twoD_Gaussian((X, Y), *coeff)
data_fitted_noise = twoD_Gaussian((X, Y), *popt)

# Calculate the extent of the plots
extent = [X.min(), X.max(), Y.min(), Y.max()]

# Plot the hist fit
plt.hist2d(x, y, bins = nbr_bins, cmap = 'terrain_r')
plt.contour(X, Y, data_fitted_hist.reshape(nbr_bins, nbr_bins), 5)
plt.scatter(0, 0, marker='x', s=100, color='r')
plt.xlim(extent[0], extent[1])  # Set the x-axis limits
plt.ylim(extent[2], extent[3])  # Set the y-axis limits
plt.gca().set_aspect('equal', adjustable='box')  # Set aspect ratio to be equal
plt.show()

# Plot the noisy data fit
fig, ax = plt.subplots(1, 1)
ax.scatter(0, 0, marker='x', s=100, color='w')
ax.imshow(data_noisy.reshape(nbr_bins, nbr_bins), cmap=plt.cm.jet, origin='lower', extent=extent)
ax.contour(X, Y, data_fitted_noise.reshape(nbr_bins, nbr_bins), 8, colors='w')
ax.set_xlim(extent[0], extent[1])  # Set the x-axis limits
ax.set_ylim(extent[2], extent[3])  # Set the y-axis limits
ax.set_aspect('equal', adjustable='box')  # Set aspect ratio to be equal
plt.show()

这是结果图：

可以清楚地看到顶部拟合（使用直方图的拟合）是不正确的，因为它没有居中（它的 sigma 参数也与另一种拟合不同）。

我想知道曲线拟合是否可能是这种拟合的错误工具？还是我做错了什么？

附注：我不知道这是否是一个线索，但在我的实际问题应用中，拟合的问题之一是它试图使高斯振幅尽可能小（如果我趋向于 0 或下限）施加界限）。

Answer 1

看来错误确实源于合并！ np.histogram2d 和 np.meshgrid 的结构似乎不一致，但这可以通过转置直方图来修复。

将 curve_fit 调用线替换为：

coeff, var_matrix = curve_fit(twoD_Gaussian,(X,Y),H.T.ravel(),p0=p0)

添加“.T”来转置直方图可以解决问题。