Generating 2D Gaussian Distributions in Python: From Independent Sampling to Multivariate Normal

Fundamental Concepts of 2D Gaussian Distribution

The 2D Gaussian distribution, also known as the bivariate normal distribution, is a crucial continuous probability distribution in probability theory and statistics. In two-dimensional space, this distribution is completely determined by a mean vector and a covariance matrix. The mean vector μ=[μ_x, μ_y] represents the center of the distribution, while the covariance matrix Σ describes the correlation between variables and their respective variances.

Independent Axis Sampling Method

When there is no correlation between the two dimensions, the covariance matrix is diagonal, and the 2D Gaussian distribution can be decomposed into the product of two independent 1D Gaussian distributions. In this case, the simplest implementation involves sampling each dimension separately:

import random

def gauss_2d_independent(mu, sigma):
    """Generate 2D Gaussian samples (independent axes)"""
    x = random.gauss(mu, sigma)
    y = random.gauss(mu, sigma)
    return (x, y)

This method is based on the mathematical principle that if X∼N(μ_x, σ_x²) and Y∼N(μ_y, σ_y²) are independent, then the joint probability density function of (X,Y) is f(x,y)=f_X(x)·f_Y(y). The random.gauss() function in the code uses either the Box-Muller transform or the ziggurat algorithm to generate standard normal random numbers, which are then transformed via linear scaling to achieve the specified mean and variance.

Multivariate Normal Distribution Method

For more general cases where correlation exists between dimensions, the multivariate normal distribution must be employed. The NumPy library provides a dedicated function for this purpose:

import numpy as np

# Define mean vector and covariance matrix
mean = np.array([0.0, 0.0])          # Shape (2,)
cov = np.array([[1.0, 0.5],          # Shape (2,2)
                [0.5, 1.0]])

# Generate 10000 samples
samples = np.random.multivariate_normal(mean, cov, 10000)

The numpy.random.multivariate_normal() function internally uses Cholesky decomposition or eigenvalue decomposition to factor the covariance matrix Σ as A·A^T, then generates samples via the linear transformation Y=μ+A·Z, where Z is a vector of standard normal random variables. This method has a time complexity of O(n³), where n is the dimensionality, but remains efficient for the 2D case.

Gaussian Kernel Matrix Generation Method

In certain image processing and signal processing applications, there is a need to directly generate 2D Gaussian kernel matrices. This approach computes Gaussian values at each grid point using an exponential function:

def make_gaussian_kernel(size, fwhm=3.0, center=None):
    """Generate a square Gaussian kernel matrix"""
    x = np.arange(0, size, 1, dtype=float)
    y = x[:, np.newaxis]
    
    if center is None:
        x0 = y0 = size // 2
    else:
        x0, y0 = center
    
    # Gaussian formula: exp(-4*ln(2)*((x-x0)^2+(y-y0)^2)/fwhm^2)
    return np.exp(-4 * np.log(2) * ((x - x0) ** 2 + (y - y0) ** 2) / fwhm ** 2)

The fwhm (full width at half maximum) parameter in the formula controls the width of the Gaussian function, related to the standard deviation σ by fwhm=2√(2ln2)·σ≈2.355·σ. This method produces a deterministic matrix of Gaussian function values rather than random samples, making it suitable for convolution kernels and filter design.

Method Comparison and Application Scenarios

Each of the three methods serves different purposes: independent axis sampling is simple and efficient for simulations with uncorrelated dimensions; multivariate normal distribution is versatile, handling arbitrary covariance structures for statistical modeling and machine learning; Gaussian kernel generation is specialized for image processing and signal filtering.

Regarding performance, for generating large sample sets, numpy.random.multivariate_normal() is highly optimized and typically faster than multiple calls to random.gauss(). Gaussian kernel generation has a computational complexity of O(n²), where n is the kernel size, making it suitable for precomputing fixed kernels.

In-Depth Mathematical Analysis

The probability density function of a 2D Gaussian distribution is:

f(x,y)=[1/(2π|Σ|^1/2)]·exp[-½·(X-μ)^TΣ^-1(X-μ)]

where X=[x,y]^T. When Σ is a diagonal matrix diag(σ_x², σ_y²), this function simplifies to the product of two 1D Gaussian functions, which is the theoretical foundation of the independent axis sampling method.

The covariance matrix Σ must be a positive semi-definite symmetric matrix, decomposable via eigenvalue decomposition as Σ=QΛQ^T, where Q is an orthogonal matrix and Λ is a diagonal matrix. Multivariate normal sampling can be achieved through the transformation Y=μ+Q√Λ·Z, where Z∼N(0,I).

Practical Application Examples

In data visualization, these methods can generate test data:

# Generate correlated data points
mean = [1.0, 2.0]
cov = [[1.0, 0.8], [0.8, 1.0]]
data = np.random.multivariate_normal(mean, cov, 500)

# Generate Gaussian filter kernel
gaussian_kernel = make_gaussian_kernel(15, fwhm=5.0)

# Apply image filtering (example)
import cv2
image = cv2.imread('input.jpg', 0)
filtered = cv2.filter2D(image, -1, gaussian_kernel)

Conclusion and Recommendations

When selecting a method for 2D Gaussian generation, consider the specific requirements: for simple independent variable sampling, two calls to random.gauss() suffice; for correlated variables, numpy.random.multivariate_normal() is recommended; for Gaussian kernels in image processing, use dedicated kernel generation functions. Understanding the mathematical principles behind each method aids in proper parameter selection and result interpretation.