Merge branch 'issue-54' into develop

* issue-54: Increment version number to 1.1.1.9015 Removed numerical computation of `qtrunc() (#54)` Adjusted unit tests (#54) Added analytical solution (#54)
ocbe-uio · Jul 3, 2024 · fc0e1f8 · fc0e1f8
2 parents 317f2b3 + 536ee4b
commit fc0e1f8
Show file tree

Hide file tree

Showing 6 changed files with 62 additions and 35 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: TruncExpFam
 Title: Truncated Exponential Family
-Version: 1.1.1.9014
+Version: 1.1.1.9015
 Date: 2024-02-26
 Authors@R:
 	c(

diff --git a/R/qtrunc.R b/R/qtrunc.R
@@ -34,29 +34,23 @@ qtrunc.generic <- function(p, ..., lower.tail, log.p) {
 }
 
 qtrunc.normal <- function(
-  p, mean = 0, sd = 1, a = -Inf, b = Inf, ..., lower.tail, log.p
+    p, mean = 0, sd = 1, a = -Inf, b = Inf, ..., lower.tail, log.p
   ) {
-  # Implements a simple bisection algorithm to find the quantile
-  lower <- rep(ifelse(a == -Inf, -.Machine$double.xmax, a), length(p))
-  upper <- rep(ifelse(b == +Inf, +.Machine$double.xmax, b), length(p))
-  q <- rowMeans(cbind(lower, upper))
-  p_q <- ptrunc.normal(q, mean, sd, a, b, lower.tail = lower.tail, log.p = log.p)
-  tol <- abs(p_q - p)
-  for (tl in seq_along(p)) {
-    iter <- 0L
-    while (tol[tl] > 1e-10 && iter < 1e9) {
-      trigger <- ifelse(lower.tail, p_q[tl] < p[tl], p_q[tl] > p[tl])
-      if (trigger) {
-        lower[tl] <- q[tl]
-        q[tl] <- mean(c(q[tl], upper[tl]))
-      } else {
-        upper[tl] <- q[tl]
-        q[tl] <- mean(c(lower[tl], q[tl]))
-      }
-      p_q[tl] <- ptrunc.normal(q[tl], mean, sd, a, b, lower.tail = lower.tail, log.p = log.p)
-      tol[tl] <- abs(p_q[tl] - p[tl])
-      iter <- iter + 1L
-    }
-  }
+  F_a <- pnorm(a, mean, sd, lower.tail, FALSE)
+  F_b <- pnorm(b, mean, sd, lower.tail, FALSE)
+  rescaled_p <- rescale_p(p, F_a, F_b, lower.tail, log.p)
+  q <- qnorm(rescaled_p, mean, sd, lower.tail, FALSE)
   return(q)
 }
+
+rescale_p <- function(p, F_a, F_b, lower.tail, log.p) {
+  if (log.p) {
+    p <- exp(p)
+  }
+  if (lower.tail) {
+    p <- p * F_b + (1 - p) * F_a
+  } else  {
+    p <- p * F_a + (1 - p) * F_b
+  }
+  return(p)
+}
diff --git a/tests/testthat/test-qtrunc-truncated-a.R b/tests/testthat/test-qtrunc-truncated-a.R
@@ -8,14 +8,22 @@ test_that("qtrunc() works as expected (normal)", {
         sg <- rchisq(1L, 5L)
         pt <- runif(i)
         if (lg) pt <- log(pt)
-        a <- qtrunc(min(pt) / 2, "normal", mean = mn, sd = sg, lower.tail = lt, log.p = lg)
-        q_trunc <- qtrunc(pt, "normal", mean = mn, sd = sg, a = a, lower.tail = lt, log.p = lg)
+        a <- qtrunc(
+          min(pt) / 2, "normal", mean = mn, sd = sg, lower.tail = lt, log.p = lg
+        )
+        q_trunc <- qtrunc(
+          pt, "normal", mean = mn, sd = sg, a = a, lower.tail = lt, log.p = lg
+        )
         q_norm <- qnorm(pt, mean = mn, sd = sg, lower.tail = lt, log.p = lg)
-        expect_length(pt, i)
         expect_length(q_trunc, i)
         for (ii in seq_along(pt)) {
           expect_gt(q_trunc[ii], q_norm[ii])
-          expect_equal(pt[ii], ptrunc(q_trunc[ii], "normal", mean = mn, sd = sg, a = a, lower.tail = lt, log.p = lg))
+          # Working back to p from q
+          ptr <- ptrunc(
+            q_trunc[ii], "normal", mean = mn, sd = sg, a = a,
+            lower.tail = lt, log.p = lg
+          )
+          expect_equal(pt[ii], ptr)
         }
       }
     }

diff --git a/tests/testthat/test-qtrunc-truncated-ab.R b/tests/testthat/test-qtrunc-truncated-ab.R
@@ -11,12 +11,20 @@ test_that("qtrunc() works as expected (normal)", {
         b <- qtrunc(max(ab), mean = mn, sd = sg, lower.tail = TRUE, log.p = FALSE)
         a <- qtrunc(min(ab), mean = mn, sd = sg, lower.tail = TRUE, log.p = FALSE)
         if (lg) pt <- log(pt)
-        q_trunc <- qtrunc(pt, "normal", mean = mn, sd = sg, a = a, b = b, lower.tail = lt, log.p = lg)
+        q_trunc <- qtrunc(
+          pt, "normal", mean = mn, sd = sg, a = a, b = b,
+          lower.tail = lt, log.p = lg
+        )
         q_norm <- qnorm(pt, mean = mn, sd = sg, lower.tail = lt, log.p = lg)
         expect_length(pt, i)
         expect_length(q_trunc, i)
         for (ii in seq_along(pt)) {
-          expect_equal(pt[ii], ptrunc(q_trunc[ii], "normal", mean = mn, sd = sg, a = a, b = b, lower.tail = lt, log.p = lg))
+          # Working back to p from q
+          ptr <- ptrunc(
+            q_trunc[ii], "normal", mean = mn, sd = sg, a = a, b = b,
+            lower.tail = lt, log.p = lg
+          )
+          expect_equal(pt[ii], ptr)
         }
       }
     }

diff --git a/tests/testthat/test-qtrunc-truncated-b.R b/tests/testthat/test-qtrunc-truncated-b.R
@@ -7,15 +7,25 @@ test_that("qtrunc() works as expected (normal)", {
         mn <- rnorm(1L, sd = 10)
         sg <- rchisq(1L, 5L)
         pt <- runif(i)
-        b <- qtrunc(max(runif(10L, pt)), "normal", mean = mn, sd = sg, lower.tail = lt, log.p = FALSE)
+        b <- qtrunc(
+          max(runif(10L, pt)), "normal", mean = mn, sd = sg,
+          lower.tail = lt, log.p = FALSE
+        )
         if (lg) pt <- log(pt)
-        q_trunc <- qtrunc(pt, "normal", mean = mn, sd = sg, b = b, lower.tail = lt, log.p = lg)
+        q_trunc <- qtrunc(
+          pt, "normal", mean = mn, sd = sg, b = b, lower.tail = lt, log.p = lg
+        )
         q_norm <- qnorm(pt, mean = mn, sd = sg, lower.tail = lt, log.p = lg)
         expect_length(pt, i)
         expect_length(q_trunc, i)
         for (ii in seq_along(pt)) {
           expect_lt(q_trunc[ii], q_norm[ii])
-          expect_equal(pt[ii], ptrunc(q_trunc[ii], "normal", mean = mn, sd = sg, b = b, lower.tail = lt, log.p = lg))
+          # Working back to p from q
+          ptr <- ptrunc(
+            q_trunc[ii], "normal", mean = mn, sd = sg, b = b,
+            lower.tail = lt, log.p = lg
+          )
+          expect_equal(pt[ii], ptr)
         }
       }
     }

diff --git a/tests/testthat/test-qtrunc-untruncated.R b/tests/testthat/test-qtrunc-untruncated.R
@@ -8,13 +8,20 @@ test_that("qtrunc() works as expected (normal)", {
         sg <- rchisq(1L, 5L)
         pt <- runif(i)
         if (lg) pt <- log(pt)
-        q_trunc <- qtrunc(pt, "normal", mean = mn, sd = sg, lower.tail = lt, log.p = lg)
+        q_trunc <- qtrunc(
+          pt, "normal", mean = mn, sd = sg, lower.tail = lt, log.p = lg
+        )
         q_norm <- qnorm(pt, mean = mn, sd = sg, lower.tail = lt, log.p = lg)
         expect_length(pt, i)
         expect_length(q_trunc, i)
         for (ii in seq_along(pt)) {
           expect_equal(q_trunc[ii], q_norm[ii])
-          expect_equal(pt[ii], ptrunc(q_trunc[ii], "normal", mean = mn, sd = sg, lower.tail = lt, log.p = lg))
+          # Working back to p from q
+          ptr <- ptrunc(
+            q_trunc[ii], "normal", mean = mn, sd = sg,
+            lower.tail = lt, log.p = lg
+          )
+          expect_equal(pt[ii], ptr)
         }
       }
     }