Density function of multivariate normal distribution for each vector x in d-dimensional space is determined as:
Where:
μ = (μ1, μ2, …, μd)T. Vector of means.
Σ. Covariance matrix with (i, j)-th element Σij = Cov(Xi, Xj). That is why Σ is symmetric positive-definite.
|Σ|. Determinant of matrix Σ.
Σ-1. Matrix inverse to Σ.
As the Σ matrix is symmetric positive-definite, it can be expanded into factors Σ = CCT (Cholesky decomposition), when the C matrix with the number of rows and columns d×d is the lower triangle. Accepting that cij is the (i, j)-th element of the C matrix, the algorithm that is used to generate required vector X with multivariate normal distribution looks as follows:
1. Execute the Cholesky decomposition for the Σ matrix to find the C matrix, which satisfies the requirement that Σ = CCT, is lower triangular and contains strictly positive diagonal entries:
Elements of the C matrix are calculated as follows:
with j < i:
2. Draw the samples Z1, Z2, …, Zd from independent and equally distributed values with the distribution N(0, 1), each sample containing n elements. Box-Muller method is used (see description of normal distribution).
3. Accept for i = 1, 2, …, d:
4. The value X = (X1, X2, …, Xd)T is returned. In matrix form X = μ + CZ.
See also: