8

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. Especially in obtaining the covariance matrix needed by numpy.random.multivariate_normal.

Franck Dernoncourt
  • 77,520
  • 72
  • 342
  • 501
Jonny
  • 1,270
  • 5
  • 19
  • 31
  • I think you are wrong thinking that ``numpy.random.multivariate_normal()`` does the same thing, because it does not give you the pdf of the distribution, it just draws random numbers from the distribution defined in the covariance matrix as well as the expectation values mu. – Nras Sep 08 '14 at 11:01
  • I see what you mean yes. Any suggestions on how to accomplish it then? – Jonny Sep 08 '14 at 11:07
  • I see your _xy_ distribution is _separable_, that is, it is the product of an _x_ Gaussian distribution times a _y_ Gaussian distribution. Perhaps that may help with Python – Luis Mendo Sep 08 '14 at 11:19

2 Answers2

8

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

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)
Richard
  • 56,349
  • 34
  • 180
  • 251
Robert Kern
  • 13,118
  • 3
  • 35
  • 32
0

I am working on a scikit called scikit-guess that contains some fast estimation routines for non-linear fits. It has a function skg.ngauss.model (also accessible as skg.ngauss_fit.model or skg.ngauss.ngauss_fit.model) which does exactly what you want. The nice thing is that it's not a PDF, so you set the amplitude out of the box:

import numpy as np
import skg.ngauss

a = 2e6
mu = 0, 0
sigma = 0.025, 0.025

x = y = np.linspace(-1, 1, 31)

cov = np.diag(sigma)**2
X = np.meshgrid(x, y)

data = skg.ngauss.model(X, a, mu, cov, axis=0)

You need to tell it axis=0 because it automatically stacks your arrays for you. To avoid passing in that argument, you could write

X = np.stack(np.meshgrid(x, y), axis=-1)

You can plot the result:

from matplotlib import pyplot as plt
plt.imshow(data)
plt.show()

enter image description here

This is not a very exciting distribution because the spread is so small that you end up with a value of ~2e-5 just one pixel away. You may want to up your sampling space to get any sort of meaningful resolution.

Note: At time of writing, the fitting function (ngauss_fit) is still buggy, but the model has been tested successfully, just not in the scikit.

Disclaimer: In case it wasn't obvious from the above, I am the author of scikit-guess.

Mad Physicist
  • 107,652
  • 25
  • 181
  • 264