可以参考R的mvtnorm包中的rmvnorm函数(生成多元正态随机数)
> rmvnorm<br />
function (n, mean = rep(0, nrow(sigma)), sigma = diag(length(mean)), <br />
method = c("svd", "chol")) <br />
{<br />
if (nrow(sigma) != ncol(sigma)) {<br />
stop("sigma must be a square matrix")<br />
}<br />
if (length(mean) != nrow(sigma)) {<br />
stop("mean and sigma have non-conforming size")<br />
}<br />
method <- match.arg(method)<br />
if (method == "svd") {<br />
ev <- eigen(sigma, sym = TRUE)$values<br />
if (!all(ev >= -sqrt(.Machine$double.eps) * abs(ev[1]))) {<br />
warning("sigma is numerically not positive definite")<br />
}<br />
sigsvd <- svd(sigma)<br />
retval <- t(sigsvd$v %*% (t(sigsvd$u) * sqrt(sigsvd$d)))<br />
}<br />
if (method == "chol") {<br />
retval <- chol(sigma, pivot = T)<br />
o <- order(attr(retval, "pivot"))<br />
retval <- retval[, o]<br />
}<br />
retval <- matrix(rnorm(n * ncol(sigma)), nrow = n) %*% retval<br />
retval <- sweep(retval, 2, mean, "+")<br />
retval<br />
}<br />
<environment: namespace:mvtnorm><br />