Add degenerate truncated normal and Student when lb = ub (fix #2)

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: TruncatedNormal
 Type: Package
 Title: Truncated Multivariate Normal and Student Distributions
-Version: 2.0
+Version: 2.0.600
 Date: 2019-06-07
 Authors@R: c(person(given="Zdravko", family="Botev", role = "aut", email = "botev@unsw.edu.au", comment = c(ORCID = "0000-0001-9054-3452")), person(given="Leo", family="Belzile", role = c("aut", "cre"), email = "belzilel@gmail.com", comment = c(ORCID = "0000-0002-9135-014X")))
 Description: A collection of functions to deal with the truncated univariate and multivariate normal and Student distributions, described in Botev (2017) <doi:10.1111/rssb.12162> and Botev and L'Ecuyer (2015) <doi:10.1109/WSC.2015.7408180>.

NEWS

@@ -1,5 +1,9 @@
 TruncatedNormal changelog
 
+version 2.1
+
+* bug fix: Cholesky decomposition with permutation now checks the arguments to ensure lb less than ub, to avoid segfault errors. In case of degenerate bounds, the conditional distribution is used with fixed components for the degenerate variables (issue #2)
+
 version 2.0
 
 * Can specify mean vector directly in the functions

R/RcppExports.R

@@ -50,7 +50,7 @@
 #' @useDynLib TruncatedNormal, .registration = TRUE
 #' @importFrom Rcpp evalCpp
 #' @export
-lnNpr <- function(a, b, check = FALSE) {
+lnNpr <- function(a, b, check = TRUE) {
     .Call(`_TruncatedNormal_lnNpr`, a, b, check)
 }
 

R/tmvnorm.R

@@ -163,10 +163,35 @@ rtmvnorm <- function(n, mu, sigma, lb, ub){
   if(missing(ub)){
     ub <- rep(Inf, d) 
   }
-  if(n == 1){
-    as.vector(mvrandn(l = lb, u = ub, Sig = sigma, n = n, mu = mu)) 
+  stopifnot(length(lb) == length(ub), length(lb) == d, lb <= ub)
+  if(!any((ub - lb) < 1e-10)){
+    if(n == 1){
+      as.vector(mvrandn(l = lb, u = ub, Sig = sigma, n = n, mu = mu)) 
+    } else{
+      t(mvrandn(l = lb, u = ub, Sig = sigma, n = n, mu = mu))
+    }
   } else{
-    t(mvrandn(l = lb, u = ub, Sig = sigma, n = n, mu = mu))
+    warning("Some variables have a degenerate distribution.")
+    ind <- which((ub - lb) >= 1e-10)
+    # check covariance matrix
+    stopifnot(isSymmetric(sigma), all(eigen(sigma, only.values = TRUE)$value > 0))
+    # compute conditional Gaussian
+      schurcomp <- function(sigma, ind) {
+        stopifnot(c(length(ind) > 0, ncol(sigma) - length(ind) > 0))
+        sigma[ind, ind, drop = FALSE] - sigma[ind, -ind, drop = FALSE] %*%
+          solve(sigma[-ind, -ind, drop = FALSE]) %*% sigma[-ind, ind, drop = FALSE]
+      }  
+      sigmap <- schurcomp(sigma, ind)
+      mup <- c(mu[ind] + sigma[ind, -ind, drop = FALSE] %*% solve(sigma[-ind, -ind, drop = FALSE]) %*% (lb[-ind] - mu[-ind]))
+      #
+      res <- matrix(0, nrow = n, ncol = d)
+      res[, -ind] <- rep(lb[-ind], each = n)
+    if(n == 1){
+      res[, ind] <- as.vector(mvrandn(l = lb[ind], u = ub[ind], Sig = sigmap, n = n, mu = mup)) 
+      return(as.vector(res))
+    } else{
+      res[, ind] <- t(mvrandn(l = lb[ind], u = ub[ind], Sig = sigmap, n = n, mu = mup))
+    }
   }
 }
 

R/tmvt.R

@@ -187,10 +187,39 @@ rtmvt <- function(n, mu, sigma, df, lb, ub){
   if(missing(ub)){
     ub <- rep(Inf, d) 
   }
-  if(n == 1){
-    as.vector(mvrandt(l = lb, u = ub, Sig = sigma, df = df, n = n, mu = mu))
+  stopifnot(length(lb) == length(ub), length(lb) == d, lb <= ub)
+  if(!any((ub - lb) < 1e-10)){
+
+    if(n == 1){
+      as.vector(mvrandt(l = lb, u = ub, Sig = sigma, df = df, n = n, mu = mu))
+    } else{
+      t(mvrandt(l = lb, u = ub, Sig = sigma, df = df, n = n, mu = mu))
+    }
   } else{
-    t(mvrandt(l = lb, u = ub, Sig = sigma, df = df, n = n, mu = mu))
+    warning("Some variables have a degenerate distribution.")
+    ind <- which((ub - lb) >= 1e-10)
+    # check covariance matrix
+    stopifnot(isSymmetric(sigma), all(eigen(sigma, only.values = TRUE)$value > 0))
+    # compute conditional Gaussian
+    schurcomp <- function(sigma, ind) {
+      stopifnot(c(length(ind) > 0, ncol(sigma) - length(ind) > 0))
+      sigma[ind, ind, drop = FALSE] - sigma[ind, -ind, drop = FALSE] %*%
+        solve(sigma[-ind, -ind, drop = FALSE]) %*% sigma[-ind, ind, drop = FALSE]
+    }  
+    sigmap <- c(df + t(lb[-ind] - mu[-ind]) %*% solve(sigma[-ind, -ind, drop = FALSE]) %*% (lb[-ind] - mu[-ind])) /
+      (df + d - length(ind)) * schurcomp(sigma, ind)
+    mup <- c(mu[ind] + sigma[ind, -ind, drop = FALSE] %*% solve(sigma[-ind, -ind, drop = FALSE]) %*% (lb[-ind] - mu[-ind]))
+    # matrix to store results
+    res <- matrix(0, nrow = n, ncol = d)
+    res[, -ind] <- rep(lb[-ind], each = n)
+    if(n == 1){
+      res[, ind] <-  as.vector(mvrandt(l = lb[ind], u = ub[ind], Sig = sigmap, 
+                                       df = df + d - length(ind), n = n, mu = mup))
+      res <- as.vector(res)
+    } else{
+      res[, ind] <- t(mvrandt(l = lb[ind], u = ub[ind], Sig = sigmap, df = df + d - length(ind), n = n, mu = mup))
+    }
+  return(res)
   }
 }