rmvnorm.Rd
Fast ways to draw multivariate normals when the variance or precision matrix is sparse.
number of observations.
mean vector.
covariance matrix (of class spam
).
precision matrix.
vector determining the mean.
the Cholesky structure of Sigma
or Q
.
arguments passed to chol
.
similar to mu
and Sigma
. Here for portability with mvtnorm::rmvnorm()
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
).
See references in chol
.
# 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