sparseMVN (version 0.2.2)

rmvn.sparse: Sample from multivariate normal distribution

Description

Efficient sampling and density calculation from a multivariate normal, when the covariance or precision matrix is sparse. These functions are designed for MVN samples of very large dimension.

Usage

rmvn.sparse(n, mu, CH, prec = TRUE)

Arguments

n

number of samples

mu

mean (numeric vector)

CH

An object of class dCHMsimpl or dCHMsuper that represents the Cholesky factorization of either the precision (default) or covariance matrix. See details.

prec

If TRUE, CH is the Cholesky decomposition of the precision matrix. If false, it is the decomposition for the covariance matrix.

Value

A matrix of samples from an MVN distribution (one in each row)

Details

This function uses sparse matrix operations to sample from a multivariate normal distribution. The user must compute the Cholesky decomposition first, using the Cholesky function in the Matrix package. This function operates on a sparse symmetric matrix, and returns an object of class dCHMsimpl or dCHMsuper (this depends on the algorithm that was used for the decomposition). This object contains information about any fill-reducing permutations that were used to preserve sparsity. The rmvn.sparse and dmvn.sparse functions use this permutation information, even if pivoting was turned off.

Examples

# NOT RUN {
   require(Matrix)
   m <- 20
   p <- 2
   k <- 4

## build sample sparse covariance matrix
   Q1 <- tril(kronecker(Matrix(seq(0.1,p,length=p*p),p,p),diag(m)))
   Q2 <- cbind(Q1,Matrix(0,m*p,k))
   Q3 <- rbind(Q2,cbind(Matrix(rnorm(k*m*p),k,m*p),Diagonal(k)))
   V <- tcrossprod(Q3)
   CH <- Cholesky(V)

   x <- rmvn.sparse(10,rep(0,p*m+k),CH, FALSE)
   y <- dmvn.sparse(x[1,],rep(0,p*m+k), CH, FALSE)

# }