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

```
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`

.

## 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
```