Draw Multivariate Normals

Description

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

Usage

rmvnorm.spam(n,mu=rep(0, nrow(Sigma)), Sigma, Rstruct=NULL, ...) 
rmvnorm.prec(n,mu=rep(0, nrow(Q)), Q, Rstruct=NULL, ...) 
rmvnorm.canonical(n, b, Q, Rstruct=NULL, ...) 

Arguments

Details

The functions rmvnorm.prec and rmvnorm.canonical do not require sparse precision matrices. For rmvnorm.spam, the differences between regular and sparse covariance matrices are too significant to be implemented here.
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).

Author(s)

Reinhard Furrer

References

See references in chol.

See Also

chol and ordering.

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")

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



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