29

I can generate Gaussian data with random.gauss(mu, sigma) function, but how can I generate 2D gaussian? Is there any function like that?

  • https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.multivariate_normal.html#scipy.stats.multivariate_normal – Moritz Jul 21 '18 at 15:55

7 Answers7

61

If you can use numpy, there is numpy.random.multivariate_normal(mean, cov[, size]).

For example, to get 10,000 2D samples:

np.random.multivariate_normal(mean, cov, 10000)

where mean.shape==(2,) and cov.shape==(2,2).

Gabriel Ziegler
  • 228
  • 2
  • 13
NPE
  • 438,426
  • 93
  • 887
  • 970
  • I am trying to draw 10000 samples from 2D distribution I created like this: data = np.random.multivariate_normal(mean,cov,(10000,10000)) but it gives memory error. Am I generating a 10000x10000 data or 2x2 data, I am confused a bit. If so, how can I draw 10000 samples from a 2D distribution? –  Oct 07 '11 at 13:34
  • I believe the correct way to get 10K 2D samples is `np.random.multivariate_normal(mean,cov,10000)`, where `mean.shape==(2,)` and `cov.shape==(2,2)`. – NPE Oct 07 '11 at 13:39
  • The above method by @NPE worked for me when I wanted to create multidimensional gaussian data. – Arnab Sinha May 26 '20 at 14:29
20

I'd like to add an approximation using exponential functions. This directly generates a 2d matrix which contains a movable, symmetric 2d gaussian.

I should note that I found this code on the scipy mailing list archives and modified it a little.

import numpy as np

def makeGaussian(size, fwhm = 3, center=None):
    """ Make a square gaussian kernel.

    size is the length of a side of the square
    fwhm is full-width-half-maximum, which
    can be thought of as an effective radius.
    """

    x = np.arange(0, size, 1, float)
    y = x[:,np.newaxis]

    if center is None:
        x0 = y0 = size // 2
    else:
        x0 = center[0]
        y0 = center[1]

    return np.exp(-4*np.log(2) * ((x-x0)**2 + (y-y0)**2) / fwhm**2)

For reference and enhancements, it is hosted as a gist here. Pull requests welcome!

giessel
  • 392
  • 3
  • 5
  • Thanks for a way to generate a matrix, that's exactly what I needed. – Phlya Oct 04 '15 at 11:27
  • I might be confused, but is the `center` actually correct? If I do `scipy.ndimage.center_of_mass(makeGaussian(10, 3, center=(4,3)))` I get `(3, 4)`. – deinonychusaur Nov 25 '15 at 11:18
  • I've also played with it a bit, the centre is indeed falsely placed – Ido_f Apr 19 '16 at 13:59
  • Check if you can simplify this using `np.outer(signal.windows.gaussian(70, 8), signal.windows.gaussian(70, 8))` like in the example on https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.fftconvolve.html – endolith Mar 11 '21 at 19:25
15

Since the standard 2D Gaussian distribution is just the product of two 1D Gaussian distribution, if there are no correlation between the two axes (i.e. the covariant matrix is diagonal), just call random.gauss twice.

def gauss_2d(mu, sigma):
    x = random.gauss(mu, sigma)
    y = random.gauss(mu, sigma)
    return (x, y)
kennytm
  • 469,458
  • 94
  • 1,022
  • 977
  • If there are no correlation between the axes, I will call random.gauss twice and I will have 2 1D gaussian dist. Then do I need to product the two 1D gaussian distribution? Or can I just put them as columns of my 2D data? Because if I need to product them together, since I have 10k data, it will cost too much. –  Oct 07 '11 at 13:46
  • @user984041: No, just treat the results as the coordinates of a 2D point. The product is the reason why this approach is valid. – kennytm Oct 07 '11 at 13:50
  • 6
    Saying to call it twice isn't a sufficient answer. You will have 2 1D arrays. A complete answer would explain how to combine two 1D arrays into a 2D array. – Octopus Jul 18 '14 at 02:14
  • 1
    @Octopus: Sampling a 2D gaussian gives you an array of 2-tuples i.e. 2×N matrix, not a 2D array (N×N matrix). I don't see how it is insufficient. – kennytm Jul 18 '14 at 07:10
  • I don't know what's meant by calling `random.gauss` twice, please modify your answer – shakram02 Jul 21 '18 at 09:05
  • 1
    @shakram02 Literally calling `random.gauss` twice . See update. – kennytm Jul 21 '18 at 14:25
6
import numpy as np

# define normalized 2D gaussian
def gaus2d(x=0, y=0, mx=0, my=0, sx=1, sy=1):
    return 1. / (2. * np.pi * sx * sy) * np.exp(-((x - mx)**2. / (2. * sx**2.) + (y - my)**2. / (2. * sy**2.)))

x = np.linspace(-5, 5)
y = np.linspace(-5, 5)
x, y = np.meshgrid(x, y) # get 2D variables instead of 1D
z = gaus2d(x, y)

Straightforward implementation and example of the 2D Gaussian function. Here sx and sy are the spreads in x and y direction, mx and my are the center coordinates.

jitter
  • 165
  • 4
  • 9
5

Numpy has a function to do this. It is documented here. Additionally to the method proposed above it allows to draw samples with arbitrary covariance.

Here is a small example, assuming ipython -pylab is started:

samples = multivariate_normal([-0.5, -0.5], [[1, 0],[0, 1]], 1000)
plot(samples[:, 0], samples[:, 1], '.')

samples = multivariate_normal([0.5, 0.5], [[0.1, 0.5],[0.5, 0.6]], 1000)
plot(samples[:, 0], samples[:, 1], '.')
Harry Cutts
  • 1,195
  • 10
  • 21
Johannes
  • 460
  • 4
  • 11
0

We can try just using the numpy method np.random.normal to generate a 2D gaussian distribution. The sample code is np.random.normal(mean, sigma, (num_samples, 2)).

A sample run by taking mean = 0 and sigma 20 is shown below :

np.random.normal(0, 20, (10,2))

>>array([[ 11.62158316,   3.30702215],
   [-18.49936277, -11.23592946],
   [ -7.54555371,  14.42238838],
   [-14.61531423,  -9.2881661 ],
   [-30.36890026,  -6.2562164 ],
   [-27.77763286, -23.56723819],
   [-18.18876597,  41.83504042],
   [-23.62068377,  21.10615509],
   [ 15.48830184, -15.42140269],
   [ 19.91510876,  26.88563983]])

Hence we got 10 samples in a 2d array with mean = 0 and sigma = 20

Abhinav Ravi
  • 713
  • 9
  • 18
0

In case someone find this thread and is looking for somethinga little more versatile (like I did), I have modified the code from @giessel. The code below will allow for asymmetry and rotation.

import numpy as np

def makeGaussian2(x_center=0, y_center=0, theta=0, sigma_x = 10, sigma_y=10, x_size=640, y_size=480):
    # x_center and y_center will be the center of the gaussian, theta will be the rotation angle
    # sigma_x and sigma_y will be the stdevs in the x and y axis before rotation
    # x_size and y_size give the size of the frame 

    theta = 2*np.pi*theta/360
    x = np.arange(0,x_size, 1, float)
    y = np.arange(0,y_size, 1, float)
    y = y[:,np.newaxis]
    sx = sigma_x
    sy = sigma_y
    x0 = x_center
    y0 = y_center

    # rotation
    a=np.cos(theta)*x -np.sin(theta)*y
    b=np.sin(theta)*x +np.cos(theta)*y
    a0=np.cos(theta)*x0 -np.sin(theta)*y0
    b0=np.sin(theta)*x0 +np.cos(theta)*y0

    return np.exp(-(((a-a0)**2)/(2*(sx**2)) + ((b-b0)**2) /(2*(sy**2))))
Messypuddle
  • 249
  • 2
  • 8