Generating 3D Gaussian distribution in Python
Question: I want to generate a Gaussian distribution in Python with the x and y
dimensions denoting position and the z dimension denoting the magnitude of a
certain quantity.

The distribution has a maximum value of 2e6 and a standard deviation
sigma=0.025.

In MATLAB I can do this with:

    
    
    x1 = linspace(-1,1,30);
    x2 = linspace(-1,1,30);
    
    mu = [0,0];
    Sigma = [.025,.025];
    
    [X1,X2] = meshgrid(x1,x2);
    F = mvnpdf([X1(:) X2(:)],mu,Sigma);
    F = 314159.153*reshape(F,length(x2),length(x1));
    surf(x1,x2,F);
    

In Python, what I have so far is:

    
    
    x = np.linspace(-1,1,30)
    y = np.linspace(-1,1,30)
    
    mu = (np.median(x),np.median(y))
    
    sigma = (.025,.025)
    

There is a Numpy function numpy.random.multivariate_normal what can supposedly
do the same as MATLAB's mvnpdf, but I am struggling to undestand the
[documentation](http://docs.scipy.org/doc/numpy/reference/generated/numpy.random.multivariate_normal.html).
Especially in obtaining the covariance matrix needed by
numpy.random.multivariate_normal.


Answer: As of scipy 0.14, you can use `scipy.stats.multivariate_normal.pdf()`:

<http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.multivariate_normal.html>

    
    
    import numpy as np
    from scipy.stats import multivariate_normal
    
    x, y = np.mgrid[-1.0:1.0:30j, -1.0:1.0:30j]
    # Need an (N, 2) array of (x, y) pairs.
    xy = np.column_stack([x.flat, y.flat])
    
    mu = np.array([0.0, 0.0])
    
    sigma = np.array([.025, .025])
    covariance = np.diag(sigma**2)
    
    z = multivariate_normal.pdf(xy, mean=mu, cov=covariance)
    
    # Reshape back to a (30, 30) grid.
    z = z.reshape(x.shape)
    

