merged in devel

.Rbuildignore

@@ -4,7 +4,7 @@
 ^\\.Rbuildignore$
 ^\\.Rhistory$
 ^\\.gitignore$
-^\\.DS\_Store
+^\\.DS_Store
 README\.md
 inst/old/*
 data-raw/*
@@ -28,3 +28,10 @@ vignettes/\_region\_.*
 inst/doc/\.install.extras
 
 
+^doc$
+^Meta$
+^_pkgdown\.yml$
+^docs$
+^pkgdown$
+
+vignettes/theory.*

.gitignore

@@ -22,3 +22,7 @@ README.html
 _region_*
 cache
 inst/doc/*
+doc
+Meta
+/doc/
+/Meta/

DESCRIPTION

@@ -2,32 +2,33 @@ Package: sparseMVN
 Type: Package
 Title: Multivariate Normal Functions for Sparse Covariance and Precision
     Matrices
-Version: 0.2.1.1
-Date: 2018-03-26
-Authors@R: person(given="Michael", family="Braun", email="[email protected]", role=c("aut","cre","cph"))
+Version: 0.2.2
+Date: 2021-10-19
+Authors@R: person(given="Michael", family="Braun", email="[email protected]", role=c("aut","cre","cph"),comment=c(ORCID="0000-0003-4774-2119"))
 Maintainer: Michael Braun <[email protected]>
-URL: http://www.smu.edu/Cox/Departments/FacultyDirectory/BraunMichael
+URL: https://braunm.github.io/sparseMVN/, https://github.com/braunm/sparseMVN/
+BugReports: https://github.com/braunm/sparseMVN/issues/
 Description: Computes multivariate normal (MVN) densities, and
     samples from MVN distributions, when the covariance or
     precision matrix is sparse.
 License: MPL (>= 2.0)
 Depends:
     R (>= 3.4.0)
 Imports:
-    Matrix (>= 1.2.12),
+    Matrix (>= 1.3),
     methods
 Suggests:
-    mvtnorm (>= 1.0.6),
-    plyr,
+    dplyr (>= 1.0),
+    tidyr (>= 1.1),
+    ggplot2 (>= 3.3),
+    forcats (>= 0.5),
+    mvtnorm (>= 1.0.6)	,
     knitr,
+    bookdown,
+    kableExtra,
     testthat,
-    dplyr (>= 0.5.0),
     scales,
-    reshape2,
-    trustOptim (>= 0.8.5),
-    xtable (>= 1.8),
-    ggplot2 (>= 2.2.1),
-    tidyr (>= 0.6.1)
-Roxygen: list(wrap=FALSE)
+    trustOptim (>= 0.8.5)
+Encoding: UTF-8
 VignetteBuilder: knitr
-RoxygenNote: 7.1.0
+RoxygenNote: 7.1.2

NEWS

@@ -1,17 +1,22 @@
+sparseMVN 0.2.2
 
+-   Added links to new pkgdown website
+-   Suggests packages used conditionally in tests
+-   Other minor fixes
 
-NEWS FILE FOR SPARSEMVN PACKAGE
+sparseMVN 0.2.1.2
 
-VERSION 0.2.1.1 (Mar. 26, 2018)
+-   Update DESCRIPTION file to be consistent with Roxygen2 >= 7.1.0
+
+sparseMVN 0.2.1.1
 
 -   Removed deprecated Matrix package functions cBind and rBind.
 
-VERSION 0.2.1
+sparseMVN 0.2.1
 
 -   New vignette.
 
-
-VERSION 0.2.0 (Feb. 5, 2015)
+sparseMVN 0.2.0
 
 -   Major overhaul of package structure, with new vignettes. Demos were
     removed in favor of vignettes.
@@ -22,8 +27,6 @@ VERSION 0.2.0 (Feb. 5, 2015)
 
 -   Now using github for development
 
-
-VERSION 0.1.0 (Nov. 1, 2013)
+sparseMVN 0.1.0
 
 -   First beta release.
-

NEWS.md

@@ -1,15 +1,23 @@
-# NEWS file for sparseMVN package
+# sparseMVN 0.2.2
 
-## VERSION 0.2.1.1 (Mar. 26, 2018)
+-  Added links to new pkgdown website
+-  Suggests packages used conditionally in tests
+-  Other minor fixes
+
+# sparseMVN 0.2.1.2
+
+-  Update DESCRIPTION file to be consistent with Roxygen2 >= 7.1.0
+
+# sparseMVN 0.2.1.1
 
 *  Removed deprecated Matrix package functions cBind and rBind.
 
 
-## VERSION 0.2.1 (May 23, 2017)
+# sparseMVN 0.2.1
 
 * New vignette.
 
-## VERSION 0.2.0 (Feb. 5, 2015)
+# sparseMVN 0.2.0
 
 * Major overhaul of package structure, with new
 vignettes.  Demos were removed in favor of vignettes.
@@ -20,7 +28,6 @@ vignettes.  Demos were removed in favor of vignettes.
 
 * Now using github for development
 
-## VERSION 0.1.0 (Nov. 1, 2013)
+# sparseMVN 0.1.0
 
 *  First beta release.
-

R/dmvn-sparse.R

@@ -0,0 +1,76 @@
+#' @rdname dmvn.sparse
+#' @title Compute density  from multivariate normal distribution
+#' @param x numeric matrix, where each row is an MVN sample.
+#' @param mu mean (numeric vector)
+#' @param CH An object of class dCHMsimpl or dCHMsuper that represents
+#' the Cholesky factorization of either the precision (default) or covariance
+#' matrix.  See details.
+#' @param prec If TRUE, CH is the Cholesky decomposition of the precision
+#' matrix.  If false, it is the decomposition for the covariance matrix.
+#' @param log If TRUE (default), returns the log density, else returns density.
+#' @return A density or log density for each row of x
+#' @section Details:
+#' This function use sparse matrix operations to 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
+#'    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)
+#'
+#' @export
+dmvn.sparse <- function(x, mu, CH, prec=TRUE, log=TRUE) {
+
+    if (is.vector(x) | (is.atomic(x) & NCOL(x)==1)) {
+        x <- matrix(x,nrow=1)
+    }
+
+    k <- length(mu)
+    n <- NROW(x)
+    if (!(k>0)) {
+        stop("mu must have positive length")
+    }
+
+    if (!(k==dim(CH)[1])) {
+        stop("dimensions of mu and CH do not conform")
+    }
+    if (k!=NCOL(x)) {
+        stop("x must have same number of columns as the length of mu")
+    }
+    if (!is.logical(prec)) {
+        stop("prec must be either TRUE or FALSE")
+    }
+    A <- expand(CH)
+    detL <- sum(log(Matrix::diag(A$L)))
+    C <- -0.918938533204672669541*k ## -k*log(2*pi)/2
+    xmu <- t(x)-mu
+    z <- as.matrix(A$P %*% xmu)
+
+    if (prec) {
+        y <- Matrix::crossprod(A$L,z)  ## L' %*% x
+        log.dens <- C + detL - Matrix::colSums(y*y)/2
+    } else {
+        y <- solve(A$L, z) ## Ly = x
+        log.dens <- C - detL - Matrix::colSums(y*y)/2
+    }
+
+    if (log) return (log.dens) else return (exp(log.dens))
+}

R/rmvn-sparse.R

@@ -1,24 +1,21 @@
 ## Copyright (C) 2013-2017 Michael Braun
-#' @title Multivariate normal functions with sparse covariance/precision matrix.
-#' @aliases dmvn.sparse rmvn.sparse
+#' @rdname rmvn.sparse
+#' @title Sample from multivariate normal distribution
+#' @aliases rmvn.sparse
 #' @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.
-#'
 #' @param n number of samples
-#' @param x numeric matrix, where each row is an MVN sample.
 #' @param mu mean (numeric vector)
 #' @param CH An object of class dCHMsimpl or dCHMsuper that represents
 #' the Cholesky factorization of either the precision (default) or covariance
 #' matrix.  See details.
 #' @param prec If TRUE, CH is the Cholesky decomposition of the precision
 #' matrix.  If false, it is the decomposition for the covariance matrix.
-#' @param log If TRUE (default), returns the log density, else returns density. 
-#'
+#' @return A matrix of samples from an MVN distribution (one in each row)
 #' @section Details:
-#' These functions uses sparse matrix operations to sample from, or compute the
-#' log density of, a multivariate normal distribution  The user must compute
+#' 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
@@ -32,14 +29,14 @@
 #'    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)
 #'
@@ -49,29 +46,22 @@ rmvn.sparse <- function(n, mu, CH, prec=TRUE) {
     if (is.na(match(class(CH),c("dCHMsimpl","dCHMsuper")))) {
         stop("CH must be an object of class 'dCHMsimpl' or 'dCHMsuper'")
     }
-
     k <- length(mu)
-
     if (!(k>0)) {
         stop("mu must have positive length")
     }
-
     if (!(n>0)) {
         stop("n must be positive")
     }
-
     if (!(k==dim(CH)[1])) {
         stop("dimensions of mu and CH do not conform")
     }
-
     if (!is.logical(prec)) {
         stop("prec must be either TRUE or FALSE")
     }
-
     x <- rnorm(n*k)
     dim(x) <- c(k,n)
-
-    A <- expand(CH)
+    A <- Matrix::expand(CH)
 
     if (prec) {
         y <- solve(Matrix::t(A$L),x) ## L'y = x
@@ -80,60 +70,7 @@ rmvn.sparse <- function(n, mu, CH, prec=TRUE) {
     }
 
     y <- as(Matrix::crossprod(A$P,y),"matrix") ## P' %*% y
-
     y <- y + mu
-
     return(t(y))
 
 }
-
-#' @rdname rmvn.sparse
-#' @export
-dmvn.sparse <- function(x, mu, CH, prec=TRUE, log=TRUE) {
-
-
-    if (is.vector(x) | (is.atomic(x) & NCOL(x)==1)) {
-        x <- matrix(x,nrow=1)
-    }
-
-    k <- length(mu)
-    n <- NROW(x)
-
-    if (!(k>0)) {
-        stop("mu must have positive length")
-    }
-
-    if (!(k==dim(CH)[1])) {
-        stop("dimensions of mu and CH do not conform")
-    }
-
-    if (k!=NCOL(x)) {
-        stop("x must have same number of columns as the length of mu")
-    }
-
-
-    if (!is.logical(prec)) {
-        stop("prec must be either TRUE or FALSE")
-    }
-
-    A <- expand(CH)
-
-    detL <- sum(log(Matrix::diag(A$L)))
-    C <- -0.918938533204672669541*k ## -k*log(2*pi)/2
-
-
-    xmu <- t(x)-mu
-
-    z <- as.matrix(A$P %*% xmu)
-
-    if (prec) {
-        y <- Matrix::crossprod(A$L,z)  ## L' %*% x
-        log.dens <- C + detL - Matrix::colSums(y*y)/2
-    } else {
-        y <- solve(A$L, z) ## Ly = x
-        log.dens <- C - detL - Matrix::colSums(y*y)/2
-    }
-
-    if (log) return (log.dens) else return (exp(log.dens))
-}
-

R/sparseMVN.R

@@ -14,5 +14,4 @@
 #' @importFrom stats rnorm
 #' @encoding UTF-8
 #' @keywords package
-NULL
-
+"_PACKAGE"