Multidimensional Fitting

[D2Hitman]

I want to fit an n-dimensional distribution with an n-dimensional gaussian.
So far i have managed to do this in 2d (see below). I am not sure how to
convert this to work in n-dimensions. Using "ravel" on the arrays is not
ideal, but optimize does not appear to work on multidimensional arrays. It
seems "meshgrid" also does not translate to nd.

import numpy as np
from scipy import optimize

#2D
#random points
N = 1000
nd = 2
a = np.random.normal(loc=[1.,2.], scale=[1.,2.], size=[N,nd])

#histogram points
bins = [10, 20]
H, e = np.histogramdd(a, bins=bins)
#find center rather than left edge
edges = []
for i in range(len(e)):
edges.append( e[:-1] + np.diff(e)/2. )

#not n-dimensional
x, y = np.meshgrid(edges[0], edges[1])
x, y = np.ravel(x), np.ravel(y)
H = np.ravel(H)

#2D gaussian
gauss2d = lambda p, x, y: p[0] * np.exp( (-0.5*(x-p[1])**2/p[2]**2) -
(0.5*(y-p[3])**2/p[4]**2) ) + p[5]
residuals = lambda p, x, y, H: H - gauss2d(p, x, y)

#Fitting
p0 = [1., 0., 1., 0., 1., 0]
plsq, cov_x, infodict, mesg, ier = optimize.leastsq(residuals, p0, args=(x,
y, H), full_output=True)

#Reading off
_x, _y = 0.091, 1.293
print '%.3f %.3f %.4f' % (_x, _y, gauss2d(plsq, _x,_y) / plsq[0])