Merge branch 'issue-54' into develop

* issue-54: Increment version number to 1.1.1.9011 Implemented `ptrunc.exp()` (#54) Implemented `ptrunc.contbern()` (#54) Implemented `ptrunc.chisq()` (#54)
ocbe-uio · Apr 9, 2024 · ad4a3b7 · ad4a3b7
2 parents 7d84b72 + 43b861f
commit ad4a3b7
Show file tree

Hide file tree

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

diff --git a/R/ptrunc.R b/R/ptrunc.R
@@ -73,13 +73,39 @@ ptrunc.poisson <- function(
   return(truncated_p(p_q, p_a, p_b, lower.tail, log.p))
 }
 
+ptrunc.chisq <- function(q, df, a = 0, b = Inf, ..., lower.tail, log.p) {
+  validate_q_a_b(q, a, b)
+  p_q <- pchisq(q, df, ncp = 0, lower.tail = TRUE, log.p)
+  p_a <- pchisq(a - 1L, df, ncp = 0, lower.tail = TRUE, log.p)
+  p_b <- pchisq(b, df, ncp = 0, lower.tail = TRUE, log.p)
+  return(truncated_p(p_q, p_a, p_b, lower.tail, log.p))
+}
+
+ptrunc.contbern <- function(q, lambda, a = 0, b = 1, ...) {
+  validate_q_a_b(q, a, b)
+  p_q <- pcontbern(q, lambda)
+  p_a <- pcontbern(a, lambda)
+  p_b <- pcontbern(b, lambda)
+  return(truncated_p(p_q, p_a, p_b, lower.tail = TRUE, log.p = FALSE))
+}
+
+ptrunc.exp <- function(q, rate = 1, a = 0, b = Inf, ..., lower.tail, log.p) {
+  validate_q_a_b(q, a, b)
+  p_q <- pexp(q, rate, lower.tail = TRUE, log.p)
+  p_a <- pexp(a, rate, lower.tail = TRUE, log.p)
+  p_b <- pexp(b, rate, lower.tail = TRUE, log.p)
+  return(truncated_p(p_q, p_a, p_b, lower.tail, log.p))
+}
+
 truncated_p <- function(p_q, p_a, p_b, lower.tail, log.p) {
+  # Handling exceptions ------------------------------------------------------
   if (!log.p && p_a == p_b) {
     return(as.numeric(lower.tail))
   }
   if (log.p && exp(p_a) == exp(p_b)) {
     return(0)
   }
+  # Usual cases --------------------------------------------------------------
   if (log.p) {
     p <- log((exp(p_q) - exp(p_a)) / (exp(p_b) - exp(p_a)))
     if (!lower.tail) {

diff --git a/tests/testthat/test-ptrunc-truncated-a.R b/tests/testthat/test-ptrunc-truncated-a.R
@@ -105,3 +105,58 @@ test_that("lower truncation works as expected (poisson)", {
     }
   }
 })
+
+test_that("lower truncation works as expected (chisq)", {
+  for (lt in c(TRUE, FALSE)) {
+    for (lg in c(FALSE, TRUE)) {
+      for (i in seq_len(10)) {
+        df <- sample(1:100, 1L)
+        a <- min(rchisq(10L, df))
+        qt <- max(rchisq(10L, df), a)
+        p_trunc <- ptrunc(
+          qt, "chisq", df, a = a, lower.tail = lt, log.p = lg
+        )
+        p_chisq <- pchisq(qt, df, ncp = 0, lower.tail = lt, log.p = lg)
+        if (!lg) {
+          expect_gte(p_trunc, 0)
+          expect_lte(p_trunc, 1)
+        } else {
+          expect_lte(p_trunc, 0)
+        }
+      }
+    }
+  }
+})
+
+test_that("lower truncation works as expected (contbern)", {
+  for (i in seq_len(10)) {
+    lambda <- runif(1L)
+    a <- runif(1L)
+    qt <- runif(1L, a, 1L)
+    p_trunc <- ptrunc(qt, "contbern", lambda, a = a)
+    p_contbern <- pcontbern(qt, lambda)
+    expect_lte(p_trunc, p_contbern)
+  }
+})
+
+test_that("lower truncation works as expected (exp)", {
+  for (lt in c(TRUE, FALSE)) {
+    for (lg in c(FALSE, TRUE)) {
+      for (i in seq_len(10)) {
+        rate <- rchisq(1L, df = 10L)
+        a <- rexp(1L, rate)
+        qt <- max(rexp(10L, rate), a)
+        p_trunc <- ptrunc(
+          qt, "exp", rate, a = a, lower.tail = lt, log.p = lg
+        )
+        p_exp <- pexp(qt, rate, lower.tail = lt, log.p = lg)
+        if (!lg) {
+          expect_gte(p_trunc, 0)
+          expect_lte(p_trunc, 1)
+        } else {
+          expect_lte(p_trunc, 0)
+        }
+      }
+    }
+  }
+})
diff --git a/tests/testthat/test-ptrunc-truncated-ab.R b/tests/testthat/test-ptrunc-truncated-ab.R
@@ -96,3 +96,61 @@ test_that("doubly-truncated ptrunc() works as expected (poisson)", {
     }
   }
 })
+
+test_that("doubly-truncated ptrunc() works as expected (chisq)", {
+  for (lt in c(TRUE, FALSE)) {
+    for (lg in c(FALSE, TRUE)) {
+      for (i in seq_len(10)) {
+        df <- sample(1:100, 1L)
+        a <- min(rchisq(10L, df))
+        b <- max(rchisq(10L, df))
+        qt <- runif(1L, a, b)
+        p_trunc <- ptrunc(
+          qt, "chisq", df, a, b, lower.tail = lt, log.p = lg
+        )
+        p_chisq <- pchisq(qt, df, ncp = 0, lower.tail = lt, log.p = lg)
+        if (!lg) {
+          expect_gte(p_trunc, 0)
+          expect_lte(p_trunc, 1)
+        } else {
+          expect_lte(p_trunc, 0)
+        }
+      }
+    }
+  }
+})
+
+test_that("doubly-truncated ptrunc() works as expected (contbern)", {
+  for (i in seq_len(10)) {
+    lambda <- runif(1L)
+    a <- runif(1L)
+    b <- runif(1L, a, 1L)
+    qt <- runif(1L, a, b)
+    p_trunc <- ptrunc(qt, "contbern", lambda, b = b)
+    p_contbern <- pcontbern(qt, lambda)
+    expect_gte(p_trunc, p_contbern)
+  }
+})
+
+test_that("doubly-truncated ptrunc() works as expected (exp)", {
+  for (lt in c(TRUE, FALSE)) {
+    for (lg in c(FALSE, TRUE)) {
+      for (i in seq_len(10)) {
+        rate <- rchisq(1L, df = 10L)
+        a <- rexp(1L, rate)
+        b <- max(rexp(10L, rate), a)
+        qt <- runif(1L, a, b)
+        p_trunc <- ptrunc(
+          qt, "exp", rate, a, b, lower.tail = lt, log.p = lg
+        )
+        p_exp <- pexp(qt, rate, lower.tail = lt, log.p = lg)
+        if (!lg) {
+          expect_gte(p_trunc, 0)
+          expect_lte(p_trunc, 1)
+        } else {
+          expect_lte(p_trunc, 0)
+        }
+      }
+    }
+  }
+})
diff --git a/tests/testthat/test-ptrunc-truncated-b.R b/tests/testthat/test-ptrunc-truncated-b.R
@@ -117,3 +117,72 @@ test_that("upper truncation works as expected (poisson)", {
     }
   }
 })
+
+test_that("upper truncation works as expected (chisq)", {
+  for (lt in c(TRUE, FALSE)) {
+    for (lg in c(FALSE, TRUE)) {
+      for (i in seq_len(10)) {
+        df <- sample(1:100, 1L)
+        b <- max(rchisq(10L, df))
+        qt <- runif(1L, 0, b)
+        p_trunc <- ptrunc(
+          qt, "chisq", df, b = b, lower.tail = lt, log.p = lg
+        )
+        p_chisq <- pchisq(qt, df, ncp = 0, lower.tail = lt, log.p = lg)
+        if (!lg) {
+          expect_gte(p_trunc, 0)
+          expect_lte(p_trunc, 1)
+          if (abs(p_trunc - p_chisq) > 1e-10) {  # adding tolerance
+            if (lt) {
+              expect_gte(p_trunc, p_chisq)
+            } else {
+              expect_lte(p_trunc, p_chisq)
+            }
+          }
+        } else {
+          expect_lte(p_trunc, 0)
+        }
+      }
+    }
+  }
+})
+
+test_that("upper truncation works as expected (contbern)", {
+  for (i in seq_len(10)) {
+    lambda <- runif(1L)
+    b <- runif(1L)
+    qt <- runif(1L, 0L, b)
+    p_trunc <- ptrunc(qt, "contbern", lambda, b = b)
+    p_contbern <- pcontbern(qt, lambda)
+    expect_gte(p_trunc, p_contbern)
+  }
+})
+
+test_that("upper truncation works as expected (exp)", {
+  for (lt in c(TRUE, FALSE)) {
+    for (lg in c(FALSE, TRUE)) {
+      for (i in seq_len(10)) {
+        rate <- rchisq(1L, df = 10L)
+        b <- rexp(1L, rate)
+        qt <- min(rexp(10L, rate), b)
+        p_trunc <- ptrunc(
+          qt, "exp", rate, b = b, lower.tail = lt, log.p = lg
+        )
+        p_exp <- pexp(qt, rate, lower.tail = lt, log.p = lg)
+        if (!lg) {
+          expect_gte(p_trunc, 0)
+          expect_lte(p_trunc, 1)
+          if (abs(p_trunc - p_exp) > 1e-10) {  # adding tolerance
+            if (lt) {
+              expect_gte(p_trunc, p_exp)
+            } else {
+              expect_lte(p_trunc, p_exp)
+            }
+          }
+        } else {
+          expect_lte(p_trunc, 0)
+        }
+      }
+    }
+  }
+})
diff --git a/tests/testthat/test-ptrunc-untruncated.R b/tests/testthat/test-ptrunc-untruncated.R
@@ -85,8 +85,60 @@ test_that("untruncated ptrunc() works as expected (poisson)", {
   }
 })
 
+test_that("untruncated ptrunc() works as expected (chisq)", {
+  for (lt in c(TRUE, FALSE)) {
+    for (lg in c(FALSE, TRUE)) {
+      for (i in seq_len(5)) {
+        df <- sample(1:10, 1L)
+        qt <- rchisq(i, df)
+        p_trunc <- ptrunc(qt, "chisq", df, lower.tail = lt, log.p = lg)
+        p_chisq <- pchisq(qt, df, lower.tail = lt, log.p = lg)
+        for (q in seq_along(qt)) {
+          if (!lg) {
+            expect_gte(p_trunc[q], 0)
+            expect_lte(p_trunc[q], 1)
+          }
+          expect_equal(p_trunc[q], p_chisq[q])
+        }
+      }
+    }
+  }
+})
+
+test_that("untruncated ptrunc() works as expected (contbern)", {
+  for (i in seq_len(5)) {
+    lambda <- runif(1)
+    qt <- rcontbern(i, lambda)
+    p_trunc <- ptrunc(qt, "contbern", lambda)
+    p_contbern <- pcontbern(qt, lambda)
+    for (q in seq_along(qt)) {
+      expect_equal(p_trunc[q], p_contbern[q])
+    }
+  }
+})
+
+test_that("untruncated ptrunc() works as expected (exp)", {
+  for (lt in c(TRUE, FALSE)) {
+    for (lg in c(FALSE, TRUE)) {
+      for (i in seq_len(5)) {
+        rate <- runif(1)
+        qt <- rexp(i, rate)
+        p_trunc <- ptrunc(qt, "exp", rate, lower.tail = lt, log.p = lg)
+        p_exp <- pexp(qt, rate, lower.tail = lt, log.p = lg)
+        for (q in seq_along(qt)) {
+          if (!lg) {
+            expect_gte(p_trunc[q], 0)
+            expect_lte(p_trunc[q], 1)
+          }
+          expect_equal(p_trunc[q], p_exp[q])
+        }
+      }
+    }
+  }
+})
+
 test_that("Basic errors are caught", {
-  for (distro in c("normal", "beta", "binomial", "poisson")) { # TODO: eventually use valid_distros
+  for (distro in c("normal", "beta", "binomial", "poisson", "chisq", "contbern", "exp")) { # TODO: eventually use valid_distros
     expect_error(ptrunc(2, distro, 1, 1, a = 3, b = 4), "must be in \\[a, b\\]")
     expect_error(ptrunc(2, distro, 1, 1, a = 0, b = 1), "must be in \\[a, b\\]")
     expect_error(ptrunc(2, distro, 1, 1, a = 3, b = 1), "a must be <= b")