Multivariate Normal Distribution: Formula & Code

Q: What is the formula for the multivariate normal distribution?

f(x) = (2π) −k/2 det(Σ) −1/2 exp[−½ (x − μ)ᵀ Σ −1 (x − μ)], valid when Σ is positive definite.

Q: What parameters define it?

The mean vector μ and the covariance matrix Σ. A k-variate normal has k means and k(k + 1)/2 distinct variance and covariance parameters.

Q: Does uncorrelated mean independent?

Only when the variables are jointly multivariate normal. Marginal normality alone is not enough.

Q: How is it different from running several univariate normal tests?

Univariate tests ignore the relationships between variables. See multivariable vs multivariate analysis and how to choose a statistical test . MANOVA, for instance, assumes multivariate normality within groups: see MANOVA in R .

The multivariate normal distribution is the joint distribution of several continuous variables in which every linear combination of them is normally distributed. It is fully described by two things: a mean vector μ and a covariance matrix Σ. This guide gives the density formula with every symbol explained, a bivariate example worked by hand, a reading table for ρ and Σ, and short Excel, R and Python recipes.

What is the multivariate normal distribution?

A random vector X with k components follows a multivariate normal (also called multivariate Gaussian) distribution when it can be written as X = μ + AZ, where Z is a vector of independent standard normal variables and A is a fixed matrix. Then Σ = AAᵀ. The univariate normal is the k = 1 case. We write X ~ N_k(μ, Σ).

The distribution matters in research because many multivariate methods (MANOVA, discriminant analysis, factor analysis, Hotelling’s T²) assume it, and because it is the standard model for correlated measurements. For the matrix itself, see our guide to the covariance matrix.

Multivariate normal distribution formula (PDF)

When Σ is positive definite, the probability density function is:

f(x) = (2π)^−k/2 · det(Σ)^−1/2 · exp[ −½ (x − μ)ᵀ Σ⁻¹ (x − μ) ]

x: the point where the density is evaluated, a vector of length k.
μ: the mean vector (length k), the location of the peak.
Σ: the k × k covariance matrix; symmetric, with variances on the diagonal and covariances off it.
det(Σ): its determinant, which scales the density down as the cloud spreads out.
Σ⁻¹: its inverse, the precision matrix.
(x − μ)ᵀ Σ⁻¹ (x − μ): the squared Mahalanobis distance D² between x and μ. See Mahalanobis distance.

This is the form used in the SciPy documentation for scipy.stats.multivariate_normal, where k is the dimension in the nonsingular case.

Worked example: a bivariate normal density

Suppose two test scores have means μ = (100, 50), standard deviations 15 and 10, and correlation ρ = 0.6. The numbers are chosen for arithmetic, not taken from a study.

Covariance: 0.6 × 15 × 10 = 90. So Σ = [[225, 90], [90, 100]].
Determinant: 225 × 100 − 90² = 22,500 − 8,100 = 14,400, so det(Σ)^1/2 = 120.
Inverse: Σ⁻¹ = (1/14,400) × [[100, −90], [−90, 225]] = [[0.006944, −0.00625], [−0.00625, 0.015625]].
Take the point x = (115, 55). Then x − μ = (15, 5).
D² = 0.006944 × 15² + 2 × (−0.00625) × 15 × 5 + 0.015625 × 5² = 1.5625 − 0.9375 + 0.3906 = 1.0156.
Exponential term: exp(−1.0156 / 2) = exp(−0.5078) = 0.6018.
Constant: 1 / (2π × 120) = 0.001326.
Density: f(115, 55) = 0.001326 × 0.6018 ≈ 0.000798.

A density is not a probability; it is a height of the surface. Probabilities come from integrating over a region. The value is checked against scipy.stats.multivariate_normal, which returns 0.000798.

Conditional distribution from the same example. Given X₁ = 115, X₂ is normal with mean 50 + (90/225)(115 − 100) = 56 and variance 100 − 90²/225 = 64 (standard deviation 8). Knowing X₁ moved the expected X₂ up by 6 and cut its variance from 100 to 64, which is the factor 1 − ρ² = 0.64.

How to interpret μ, Σ and ρ

Feature of Σ	What the contours look like
ρ = 0, equal variances	Circles centred on μ
ρ = 0, unequal variances	Ellipses with axes parallel to the variable axes
ρ > 0	Ellipses tilted so both variables rise together
ρ < 0	Ellipses tilted so one rises as the other falls
\|ρ\| near 1	Very thin ellipses; Σ approaches singular (det(Σ) near 0)

The contours are the sets of points with equal D². In k dimensions, D² of a normal vector follows a chi-square distribution with k degrees of freedom, which is why points with D² beyond the chi-square 95th percentile (5.99 for k = 2) lie outside the central 95% ellipse. The worked point has D² = 1.02, inside it (p = 0.60).

Properties of the multivariate normal distribution

Marginals and linear transformations

Every subset of components is itself multivariate normal, and any linear transformation AX + b of a multivariate normal vector is again multivariate normal with mean Aμ + b and covariance AΣAᵀ.

Zero correlation and independence

For jointly normal variables, zero covariance means independence. This does not hold in general: two variables that are each normal but not jointly normal can be uncorrelated and still dependent. Marginal normality of each variable does not prove joint normality.

Conditional distributions

Conditioning a multivariate normal on some components gives another multivariate normal. The conditional mean is linear in the conditioning values, and the conditional covariance does not depend on them, which is why linear regression and the multivariate normal are so closely tied.

Assumptions and common mistakes

Checking each variable separately. Univariate normality (for example a Shapiro-Wilk test) is necessary but not sufficient for joint normality.
Ignoring the joint check. Common options are Mardia’s tests of multivariate skewness and kurtosis (Mardia, 1970, Biometrika 57(3), 519–530) and a chi-square Q-Q plot of the squared Mahalanobis distances against chi-square quantiles with k degrees of freedom.
Treating a density as a probability. Density values can exceed 1 when variances are small.
Using a Σ that is not valid. A covariance matrix must be symmetric and positive semidefinite; hand-typed matrices often are not.
Over-trusting tests in large samples. Formal normality tests reject trivial departures when n is large, so pair them with plots.

Maximum likelihood estimation (MLE)

For a sample of n independent observations, the maximum likelihood estimate of μ is the sample mean vector. The MLE of Σ is the sample covariance with divisor n, not n − 1: Σ̂ = (1/n) Σᵢ (xᵢ − x̄)(xᵢ − x̄)ᵀ. Most software reports the unbiased version with divisor n − 1, so the two differ by the factor n/(n − 1).

Singular multivariate normal distribution

If Σ is only positive semidefinite (det Σ = 0, for example when one variable is an exact linear function of another), no density exists in the full k-dimensional space; the distribution lives on a lower-dimensional subspace. The SciPy documentation handles this with allow_singular=True, using the pseudo-determinant and pseudo-inverse, and defines k there as the rank of Σ. Otherwise the covariance must be strictly positive definite.

How to report it (APA example)

Report the parameters and how you checked the assumption. A template, using the example’s numbers:

The two scores were modeled as bivariate normal with mean vector μ = (100, 50) and covariance matrix Σ = [[225, 90], [90, 100]] (r = .60). Multivariate normality was assessed with Mardia’s skewness and kurtosis tests and a chi-square Q-Q plot of squared Mahalanobis distances.

Italicize statistical symbols such as r and p, and drop the leading zero for values that cannot exceed 1. Replace the numbers and the test outcome with your own results.

How to calculate it in Excel, R or Python

Excel

There is no built-in multivariate normal function. With the deviation row (x − μ) in A5:B5 and Σ in A1:B2, an array formula gives the density:

=EXP(-0.5*MMULT(MMULT(A5:B5,MINVERSE(A1:B2)),TRANSPOSE(A5:B5)))/(2*PI()*SQRT(MDETERM(A1:B2)))

R

library(mvtnorm)
S <- matrix(c(225, 90, 90, 100), 2)
dmvnorm(c(115, 55), mean = c(100, 50), sigma = S)   # density
rmvnorm(1000, mean = c(100, 50), sigma = S)         # simulate

Python (SciPy and NumPy)

import numpy as np
from scipy.stats import multivariate_normal as mvn
S = [[225, 90], [90, 100]]
mvn(mean=[100, 50], cov=S).pdf([115, 55])   # 0.000798
rng = np.random.default_rng(1)
X = rng.multivariate_normal([100, 50], S, size=1000)
np.cov(X, rowvar=False, bias=True)          # MLE of Sigma

Frequently asked questions

What is the formula for the multivariate normal distribution?

f(x) = (2π)^−k/2 det(Σ)^−1/2 exp[−½ (x − μ)ᵀ Σ⁻¹ (x − μ)], valid when Σ is positive definite.

What parameters define it?