Multivariate normal functions with sparse covariance/precision matrix.

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

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

dmvn.sparse(x, 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.

x

numeric matrix, where each row is an MVN sample.

Details

These functions uses sparse matrix operations to sample from, or compute the log density of, 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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
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)
 ##  print(x)

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

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.