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

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
– Nras, Commented Sep 8, 2014 at 11:01
I see what you mean yes. Any suggestions on how to accomplish it then? — Jonny
– Jonny, Commented Sep 8, 2014 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
– Luis Mendo, Commented Sep 8, 2014 at 11:19

Richard · Accepted Answer · 2020-05-20 17:26:52Z

9

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)

edited May 20, 2020 at 17:26

Richard

62.7k40 gold badges198 silver badges277 bronze badges

answered Sep 8, 2014 at 11:25

Robert Kern

13.5k3 gold badges37 silver badges32 bronze badges

Sign up to request clarification or add additional context in comments.

Comments

Mad Physicist · Accepted Answer · 2021-02-15 03:06:56Z

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()

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.

Collectives™ on Stack Overflow

Generating 3D Gaussian distribution in Python

2 Answers 2

Comments

Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related