The Gaussian (normal) distribution is the workhorse of machine learning. It’s the default model for noise, the likelihood and prior in linear regression, the components of a mixture model, and the limiting shape of averages (the central limit theorem). Its algebra is extraordinarily convenient: marginals, conditionals, products, and linear transforms of Gaussians are all Gaussian, giving closed-form answers where other distributions need approximation.

A real-life example: measurement noise

Weigh the same object 100 times on a digital scale and the readings scatter around the true weight in a symmetric bell shape — a Gaussian. This isn’t a coincidence: the central limit theorem says that a sum of many small independent errors (air currents, vibration, rounding) tends toward a Gaussian regardless of their individual shapes. That’s why “assume Gaussian noise” is the sensible default across science and ML.

The univariate Gaussian

A single random variable is Gaussian, $X \sim \mathcal{N}(\mu, \sigma^2)$, with density

The mean $\mu$ centers the bell; the variance $\sigma^2$ sets its width. About 68% of mass lies within $1\sigma$, 95% within $2\sigma$.

The multivariate Gaussian

In $D$ dimensions, $\mathbf{x} \sim \mathcal{N}(\boldsymbol\mu, \Sigma)$ with mean vector $\boldsymbol\mu$ and covariance matrix $\Sigma$:

The covariance $\Sigma$ shapes the contours into ellipses: diagonal $\Sigma$ gives axis-aligned ellipses; off-diagonal (correlation) tilts them. The special case $\boldsymbol\mu = \mathbf{0}$, $\Sigma = I$ is the standard normal.

Watch covariance shape the Gaussian

The contours of a 2-D Gaussian are ellipses. As the correlation in $\Sigma$ sweeps, the ellipse tilts — positive correlation stretches it along the diagonal, zero makes it a circle, negative tilts it the other way:

Why Gaussians are so convenient: closure

The Gaussian’s superpower is closure — Gaussian operations return Gaussians:

• Marginals of a joint Gaussian are Gaussian (just read off the relevant $\boldsymbol\mu$, $\Sigma$ blocks).
• Conditionals are Gaussian: $p(\mathbf{x}\mid\mathbf{y}) = \mathcal{N}(\boldsymbol\mu_{x\mid y}, \Sigma_{x\mid y})$ with closed-form mean/covariance — the core of Kalman filters and Gaussian processes.
• Products of Gaussian densities are (scaled) Gaussian — used to multiply prior × likelihood.
• Sums of independent Gaussians are Gaussian: means and covariances add.
• Linear/affine transforms stay Gaussian: if $\mathbf{x} \sim \mathcal{N}(\boldsymbol\mu, \Sigma)$ then $A\mathbf{x} + \mathbf{b} \sim \mathcal{N}(A\boldsymbol\mu + \mathbf{b}, A\Sigma A^\top)$.

  flowchart TD
  G["Gaussian N(μ, Σ)"] --> M["marginal → Gaussian"]
  G --> C["conditional → Gaussian"]
  G --> P["product of densities → Gaussian"]
  G --> S["sum of independents → Gaussian"]
  G --> L["linear transform Ax+b → N(Aμ+b, AΣAᵀ)"]
  L -.-> SAMP["sampling: μ + L z, Σ = LLᵀ (Cholesky)"]

That last property gives the sampling recipe: to draw from $\mathcal{N}(\boldsymbol\mu, \Sigma)$, take $\mathbf{z} \sim \mathcal{N}(\mathbf{0}, I)$ and return $\boldsymbol\mu + L\mathbf{z}$ where $\Sigma = LL^\top$ (the Cholesky factor).

NumPy

import numpy as np
rng = np.random.default_rng(0)
 
# univariate density
def normal_pdf(x, mu, sigma):
    return np.exp(-0.5*((x-mu)/sigma)**2) / (sigma*np.sqrt(2*np.pi))
print("N(0,1) at x=0:", round(normal_pdf(0, 0, 1), 4))   # 0.3989
 
# sample a 2-D Gaussian via Cholesky:  mu + L z
mu = np.array([1.0, -2.0])
Sigma = np.array([[2.0, 1.2],
                  [1.2, 1.5]])
L = np.linalg.cholesky(Sigma)
z = rng.standard_normal((2, 100000))
samples = mu[:, None] + L @ z
print("empirical mean:", np.round(samples.mean(axis=1), 2))
print("empirical cov:\n", np.round(np.cov(samples), 2))
 
# linear transform stays Gaussian: y = A x + b
A = np.array([[0.5, 0.0], [0.0, 2.0]])
print("transformed cov = A Σ Aᵀ:\n", np.round(A @ Sigma @ A.T, 2))

import numpy as np
rng = np.random.default_rng(0)
 
# univariate density
def normal_pdf(x, mu, sigma):
    return np.exp(-0.5*((x-mu)/sigma)**2) / (sigma*np.sqrt(2*np.pi))
print("N(0,1) at x=0:", round(normal_pdf(0, 0, 1), 4))   # 0.3989
 
# sample a 2-D Gaussian via Cholesky:  mu + L z
mu = np.array([1.0, -2.0])
Sigma = np.array([[2.0, 1.2],
                  [1.2, 1.5]])
L = np.linalg.cholesky(Sigma)
z = rng.standard_normal((2, 100000))
samples = mu[:, None] + L @ z
print("empirical mean:", np.round(samples.mean(axis=1), 2))
print("empirical cov:\n", np.round(np.cov(samples), 2))
 
# linear transform stays Gaussian: y = A x + b
A = np.array([[0.5, 0.0], [0.0, 2.0]])
print("transformed cov = A Σ Aᵀ:\n", np.round(A @ Sigma @ A.T, 2))

N(0,1) at x=0: 0.3989
empirical mean: [ 1.  -2. ]
empirical cov:
 [[2.   1.19]
 [1.19 1.5 ]]
transformed cov = A Σ Aᵀ:
 [[0.5 1.2]
 [1.2 6. ]]

N(0,1) at x=0: 0.3989
empirical mean: [ 1.  -2. ]
empirical cov:
 [[2.   1.19]
 [1.19 1.5 ]]
transformed cov = A Σ Aᵀ:
 [[0.5 1.2]
 [1.2 6. ]]

Why this matters for ML

• Linear regression (Chapter 9) uses a Gaussian likelihood and prior; closure gives a closed-form Gaussian posterior.
• Gaussian processes, Kalman filters, variational autoencoders, PPCA all exploit Gaussian marginals/conditionals for tractable inference.
• Gaussian noise is the default assumption; the CLT justifies it whenever many small effects combine.

🧪 Try It Yourself

Exercise 1 – The Gaussian density

Exercise 2 – Sample via Cholesky

Exercise 3 – Linear transform stays Gaussian

Recap

• The Gaussian $\mathcal{N}(\mu, \sigma^2)$ / $\mathcal{N}(\boldsymbol\mu, \Sigma)$ is the bell curve; covariance shapes its elliptical contours.
• It is closed under marginalization, conditioning, products, sums, and linear transforms — all return Gaussians with closed-form parameters.
• Sample via $\boldsymbol\mu + L\mathbf{z}$ (Cholesky); the CLT explains its ubiquity.

Next: distributions that pair nicely for Bayesian updates — Conjugacy and the Exponential Family.