Steps to sample from a multivariate Gaussian (Normal) distribution with Python code
A widely used method for drawing (sampling) a random vector x from the N-dimensional multivariate normal distribution with mean vector μ and covariance matrixΣ works as follows:
Find any real matrix A such that AAT = Σ. When Σ is positive-definite, the Cholesky decomposition is typically used, and the extended form of this decomposition can always be used (as the covariance matrix may be only positive semi-definite) in both cases a suitable matrix A is obtained. An alternative is to use the matrix A = UΛ½ obtained from a spectral decompositionΣ = UΛU−1 of Σ. The former approach is more computationally straightforward but the matrices A change for different orderings of the elements of the random vector, while the latter approach gives matrices that are related by simple re-orderings. In theory both approaches give equally good ways of determining a suitable matrix A, but there are differences in computation time.