Draw Multivariate Normals

Fast ways to draw multivariate normals when the variance or precision matrix is sparse.

Usage

rmvnorm(n, mu=rep.int(0, dim(Sigma)[1]), Sigma, ..., mean, sigma)
rmvnorm.spam(n, mu=rep.int(0, dim(Sigma)[1]), Sigma, Rstruct=NULL, ..., mean, sigma)
rmvnorm.prec(n, mu=rep.int(0, dim(Q)[1]), Q, Rstruct=NULL, ...)
rmvnorm.canonical(n, b, Q, Rstruct=NULL, ...)

Arguments

n: number of observations.
mu: mean vector.
Sigma: covariance matrix (of class spam).
Q: precision matrix.
b: vector determining the mean.
Rstruct: the Cholesky structure of Sigma or Q.
...: arguments passed to chol.
mean,sigma: similar to mu and Sigma. Here for portability with mvtnorm::rmvnorm()

Details

All functions rely on a Cholesky factorization of the covariance or precision matrix.
The functions rmvnorm.prec and rmvnorm.canonical do not require sparse precision matrices Depending on the the covariance matrix Sigma, rmvnorm or rmvnorm.spam is used. If wrongly specified, dispatching to the other function is done.
Default mean is zero. Side note: mean is added via sweep() and no gain is accieved by distinguishing this case.
Often (e.g., in a Gibbs sampler setting), the sparsity structure of the covariance/precision does not change. In such setting, the Cholesky factor can be passed via Rstruct in which only updates are performed (i.e., update.spam.chol.NgPeyton instead of a full chol).

References

See references in chol.

Author

Reinhard Furrer

Examples

# Generate multivariate from a covariance inverse:
# (usefull for GRMF)
set.seed(13)
n <- 25    # dimension
N <- 1000  # sample size
Sigmainv <- .25^abs(outer(1:n,1:n,"-"))
Sigmainv <- as.spam( Sigmainv, eps=1e-4)


Sigma <- solve( Sigmainv)  # for verification
iidsample <- array(rnorm(N*n),c(n,N))

mvsample <- backsolve( chol(Sigmainv), iidsample)
norm( var(t(mvsample)) - Sigma, type="m")
#> [1] 0.1326448

# compare with:
mvsample <- backsolve( chol(as.matrix( Sigmainv)), iidsample, n)
   #### ,n as patch
norm( var(t(mvsample)) - Sigma, type="m")
#> [1] 0.1326447


# 'solve' step by step:
b <- rnorm( n)
R <- chol(Sigmainv)
norm( backsolve( R, forwardsolve( R, b))-
      solve( Sigmainv, b) )
#> [1] 0
norm( backsolve( R, forwardsolve( R, diag(n)))- Sigma )
#> [1] 0