Fix error in sanity checks

lbelzile · May 16, 2020 · 01ada6c · 01ada6c
1 parent 2029646
commit 01ada6c
Show file tree

Hide file tree

Showing 3 changed files with 16 additions and 20 deletions.
diff --git a/R/tmvnorm.R b/R/tmvnorm.R
@@ -157,7 +157,7 @@ ptmvnorm <- function(q, mu, sigma, lb, ub, log = FALSE, type = c("mc", "qmc"), B
       pb <- switch(type,
                         mc = mvNcdf(l = lb - mu, u = pmin(ub, q[i,]) - mu, Sig = sigma, n = B)$prob,
                         qmc = mvNqmc(l = lb - mu, u = pmin(ub, q[i,]) - mu, Sig = sigma, n = B)$prob)
-      prob[i] <- ifelse(log, pmin(0, log(prob[i]) - log(kst)), pmin(1, pmax(0, prob[i]/kst)))
+      prob[i] <- ifelse(log, pmin(0, log(pb) - log(kst)), pmin(1, pb/kst))
     }
   }
   return(prob)

diff --git a/R/tmvt.R b/R/tmvt.R
@@ -140,26 +140,22 @@ ptmvt <- function(q, mu, sigma, df, lb, ub, type = c("mc", "qmc"), log = FALSE,
   }
   stopifnot(all(lb < ub))
   prob <- rep(0, nrow(q))
-
+  kst <- switch(type,
+                mc = mvTcdf(l = lb - mu, u = ub - mu, df = df, Sig = sigma, n = B)$prob,
+                qmc = mvTqmc(l = lb - mu, u = ub - mu, df = df, Sig = sigma, n = B)$prob)
   for(i in 1:nrow(q)){
     if(all(q[i,] >= ub)){
-      prob[i] <- 1
+      prob[i] <- ifelse(log, 0, 1)
     } else if(any(q[i,] <= lb)){
-      prob[i] <- 0
+      prob[i] <- ifelse(log, -Inf, 0)
     } else{
-      prob[i] <- switch(type,
-                        mc = mvTcdf(l = lb - mu, u = pmin(ub, q[i,]) - mu, df = df, Sig = sigma, n = B)$prob,
-                        qmc = mvTqmc(l = lb - mu, u = pmin(ub, q[i,]) - mu, df = df, Sig = sigma, n = B)$prob)
+      pb <- switch(type,
+                   mc = mvTcdf(l = lb - mu, u = pmin(ub, q[i,]) - mu, df = df, Sig = sigma, n = B)$prob,
+                   qmc = mvTqmc(l = lb - mu, u = pmin(ub, q[i,]) - mu, df = df, Sig = sigma, n = B)$prob)
+      prob[i] <- ifelse(log, pmin(0, log(pb) - log(kst)), pmin(1, pb/kst))
     }
   }
-  kst <- switch(type,
-                mc = mvTcdf(l = lb - mu, u = ub - mu, df = df, Sig = sigma, n = B)$prob,
-                qmc = mvTqmc(l = lb - mu, u = ub - mu, df = df, Sig = sigma, n = B)$prob)
-  if(log){
-    return(pmin(0, log(prob) - log(kst)))
-  } else{
-    return(pmin(1,pmax(0, prob/kst)))
-  }
+  return(prob)
 }
 
 #' Random number generator for the truncated multivariate Student distribution.

diff --git a/tests/testthat/test-wrappers.R b/tests/testthat/test-wrappers.R
@@ -8,7 +8,7 @@ df2 <- 300
 x = c(100, 50)
 
 testthat::test_that("Truncated normal DF (MC versus QMC) give similar answers", {
-  skip_on_cran()
+  testthat::skip_on_cran()
   testthat::expect_equal(ptmvnorm(x, mu=mu, sigma=sigma, lb=lb, ub=ub, log=FALSE, type="qmc"),
                          ptmvnorm(x, mu=mu, sigma=sigma, lb=lb, ub=ub, log=FALSE, type="mc"), 
                          tolerance = 1e-4)
@@ -18,7 +18,7 @@ testthat::test_that("Truncated normal DF (MC versus QMC) give similar answers",
 })
 
 testthat::test_that("Student with large df gives same answer as normal", {
-  skip_on_cran()
+  testthat::skip_on_cran()
   testthat::expect_equal(ptmvnorm(x, B = 1e6, sigma=sigma, lb=lb, ub=ub, log=FALSE, type="qmc"),
                          ptmvt(x, B = 1e6, sigma=sigma, lb=lb, ub=ub, df = 300, log=FALSE, type="qmc"), 
                          tolerance = 2e-3)
@@ -55,7 +55,7 @@ D <- 10
 muV <- 1:D
 Smat <- diag(0.5, D) + matrix(0.5, D, D)
 testthat::test_that("Expectation of (truncated) elliptical distributions", {
-  skip_on_cran()
+  testthat::skip_on_cran()
   testthat::expect_equal(colMeans(rtmvnorm(n = B, sigma = Smat)),
                          rep(0, D), tolerance = 5/sqrt(B))
   testthat::expect_equal(colMeans(rtmvnorm(n = B, mu = muV, sigma = Smat)),
@@ -72,7 +72,7 @@ testthat::test_that("Expectation of (truncated) elliptical distributions", {
 lb <- rnorm(n = D, mean = 0, sd = 10)
 ub <- lb + rgamma(D, shape = 4, rate = 1)
 testthat::test_that("Bounds of simulated variables", {
-  skip_on_cran()
+  testthat::skip_on_cran()
   testthat::expect_true(isTRUE(all(apply(rtmvnorm(n = 1e4, lb = lb, ub = ub, mu = muV, 100*Smat), 2, min) > lb)))
   testthat::expect_true(isTRUE(all(apply(rtmvnorm(n = 1e4, lb = lb, ub = ub, mu = muV, 100*Smat), 2, min) < ub)))
   testthat::expect_true(isTRUE(all(apply(rtmvt(n = 1e4, df = 3, lb = lb, ub = ub, mu = muV, 100*Smat), 2, min) > lb)))
@@ -107,7 +107,7 @@ testthat::test_that("Untruncated density agrees with that in the mvtnorm package
 d <- 15
 sigma <- 0.5 * (diag(d) + matrix(1, d, d))
 testthat::test_that("Known probability", {
-  skip_on_cran()
+  testthat::skip_on_cran()
   testthat::expect_equivalent((d+1)*pmvnorm(sigma = sigma, lb = rep(0, d), type = "qmc", B = B),
   1, tolerance = 1/sqrt(B))
   })