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.

1 2 3 | ```
rmvn.sparse(n, mu, CH, prec = TRUE)
dmvn.sparse(x, mu, CH, prec = TRUE)
``` |

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

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.

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

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

All documentation is copyright its authors; we didn't write any of that.