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