如何将二维轮廓图拟合到二维高斯分布

Question

我有一个二维轮廓图，我想用二维高斯函数来拟合它。这是我用来绘制二维轮廓图的脚本。

import numpy as np
from pylab import *
from scipy.stats import kde
x = np.genfromtxt("deltaDEC.dat",delimiter="\n")
y = np.genfromtxt("cosDEC.dat",delimiter="\n")
n = len(x)
H, xedges, yedges = np.histogram2d(x, y, range=[[-40,40], [-40,40]], bins=(50, 50))
extent = [yedges[0], yedges[-1], xedges[0], xedges[-1]]
levels = (400, 200, 100, 50, 20)
cset = contour(H, levels, origin='lower',colors=['black', 'pink','green','blue','red'],linewidths=(1.9, 1.6, 1.5, 1.4),extent=extent)
clabel(cset, inline=1, fontsize=10, fmt='%1.0i')
ylim(-10, 10)
xlim(-40, 40)
xlabel('delta_RA/cos(DEC)')
ylabel('delta_DEC')
for c in cset.collections:
  c.set_linestyle('solid')
colorbar()
show()

那我该怎么拟合呢？

编辑

就像这个脚本，它是用来将高斯函数拟合到一维直方图的，但我想把高斯函数拟合到二维直方图上。

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> import matplotlib.mlab as mlab
>>> from scipy.stats import norm
>>> data = np.loadtxt('delta DEC".txt')
>>> (mu,sigma) = norm.fit(data)
>>> plt.figure(1)
<matplotlib.figure.Figure object at 0x26bc350>
>>> n, bins, patches=plt.hist(data, 100, normed=True, histtype='step', facecolor='green')
>>> y = mlab.normpdf(bins, mu, sigma)
>>> plt.plot(bins, y, 'r--', linewidth=2)
[<matplotlib.lines.Line2D object at 0x31d0a10>]
>>> plt.xlabel('DEC difference[arcsec]')
<matplotlib.text.Text object at 0x31cb910>
>>> plt.ylabel('Probability')
<matplotlib.text.Text object at 0x2e8c3d0>
>>> plt.title('Gaussian distribution')
<matplotlib.text.Text object at 0x31d66d0>
>>> plt.grid(True)
>>> plt.show()

Answer 1

如果你想把一个二维的高斯分布（也就是钟形曲线）适配到你的数据 x 和 y 上，其实非常简单（Numpy里已经内置了计算均值和协方差的最大似然估计）——比如说，下面是针对随机生成的 x 和 y 数据的例子。

>>> import numpy as np
>>> x, y = np.random.randn(2, 100)
>>> data = np.vstack([x, y])

>>> np.mean(data, axis=1)
array([ 0.01154114, -0.01544327])

>>> np.cov(data)
array([[ 1.19047626, -0.11507689],
       [-0.11507689,  0.95112915]])