Skip to content

Commit

Permalink
Merge branch 'issue-54' into develop
Browse files Browse the repository at this point in the history 
* issue-54:
  Added tests, removed unused validation (#54)
  Increment version number to 1.1.1.9012
  Fixed unit tests (#54)
  Standardized test vector length (#54)
  Added tests for `p_trunc` length + fixes (#54)
  Added tests for `qt` length (#54)
  Implemented `ptrunc.invgauss()` (#54)
  Implemented `ptrunc.invgamma()` (#54)
  Implemented `ptrunc.gamma` (#54)
  • Loading branch information
@wleoncio
wleoncio committed Apr 16, 2024
2 parents ad4a3b7 + e873054 commit 5bb4d85
Show file tree
Hide file tree
Showing 6 changed files with 676 additions and 192 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Package: TruncExpFam
Title: Truncated Exponential Family
Version: 1.1.1.9011
Version: 1.1.1.9012
Date: 2024-02-26
Authors@R:
c(
Expand Down
39 changes: 33 additions & 6 deletions R/ptrunc.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,14 +97,41 @@ ptrunc.exp <- function(q, rate = 1, a = 0, b = Inf, ..., lower.tail, 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))
ptrunc.gamma <- function(
q, shape, rate = 1, scale = 1 / rate, a = 0, b = Inf, ..., lower.tail, log.p
) {
validate_q_a_b(q, a, b)
if (!missing(rate) && !missing(scale)) {
stop("specify 'rate' or 'scale' but not both")
}
if (log.p && exp(p_a) == exp(p_b)) {
return(0)
p_q <- pgamma(q, shape, scale = scale, lower.tail = TRUE, log.p = log.p)
p_a <- pgamma(a, shape, scale = scale, lower.tail = TRUE, log.p = log.p)
p_b <- pgamma(b, shape, scale = scale, lower.tail = TRUE, log.p = log.p)
return(truncated_p(p_q, p_a, p_b, lower.tail, log.p))
}

ptrunc.invgamma <- function(
q, shape, rate = 1, scale = 1 / rate, a = 0, b = Inf, ..., lower.tail, log.p
) {
validate_q_a_b(q, a, b)
if (!missing(rate) && !missing(scale)) {
stop("specify 'rate' or 'scale' but not both")
}
p_q <- pinvgamma(q, shape, scale = scale, lower.tail = TRUE, log.p = log.p)
p_a <- pinvgamma(a, shape, scale = scale, lower.tail = TRUE, log.p = log.p)
p_b <- pinvgamma(b, shape, scale = scale, lower.tail = TRUE, log.p = log.p)
return(truncated_p(p_q, p_a, p_b, lower.tail, log.p))
}

ptrunc.invgauss <- function(q, m, s, a = 0, b = Inf, ...) {
validate_q_a_b(q, a, b)
p_q <- pinvgauss(q, m, s)
p_a <- ifelse(a == 0, 0, pinvgauss(a, m, s))
p_b <- ifelse(b == Inf, 1, pinvgauss(b, m, s))
return(truncated_p(p_q, p_a, p_b, lower.tail = TRUE, log.p = FALSE))
}

truncated_p <- function(p_q, p_a, p_b, lower.tail, log.p) {
# Usual cases --------------------------------------------------------------
if (log.p) {
p <- log((exp(p_q) - exp(p_a)) / (exp(p_b) - exp(p_a)))
Expand Down
233 changes: 173 additions & 60 deletions tests/testthat/test-ptrunc-truncated-a.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,25 +5,29 @@ test_that("lower truncation works as expected (normal)", {
lg <- FALSE
for (lt in c(TRUE, FALSE)) {
for (lg in c(FALSE, TRUE)) {
for (i in seq_len(1)) {
for (i in seq_len(5)) {
mn <- rnorm(1L, sd = 10)
sg <- rchisq(1L, 5L)
qt <- rnorm(1L, mn, sg)
a <- qt - rchisq(1L, 5L)
qt <- rnorm(i, mn, sg)
a <- min(qt) - rchisq(1L, 5L)
p_trunc <- ptrunc(
qt, lower.tail = lt, log.p = lg, mean = mn, sd = sg, a = a
)
p_norm <- pnorm(qt, lower.tail = lt, log.p = lg, mean = mn, sd = sg)
if (!lg) {
expect_gte(p_trunc, 0)
expect_lte(p_trunc, 1)
if (lt) {
expect_lte(p_trunc, p_norm)
expect_length(qt, i)
expect_length(p_trunc, i)
for (q in seq_along(qt)) {
if (!lg) {
expect_gte(p_trunc[q], 0)
expect_lte(p_trunc[q], 1)
if (lt) {
expect_lte(p_trunc[q], p_norm[q])
} else {
expect_gte(p_trunc[q], p_norm[q])
}
} else {
expect_gte(p_trunc, p_norm)
expect_lte(p_trunc[q], 0)
}
} else {
expect_lte(p_trunc, 0)
}
}
}
Expand All @@ -33,27 +37,29 @@ test_that("lower truncation works as expected (normal)", {
test_that("lower truncation works as expected (beta)", {
for (lt in c(TRUE, FALSE)) {
for (lg in c(FALSE, TRUE)) {
for (i in seq_len(10)) {
for (i in seq_len(5)) {
shp1 <- sample(1:10, 1L)
shp2 <- sample(1:10, 1L)
list2env(
setNames(as.list(sort(rbeta(2L, shp1, shp2))), c("a", "qt")),
envir = .GlobalEnv
)
a <- rbeta(1L, shp1, shp2)
qt <- replicate(i, max(rbeta(10L, shp1, shp2), a))
p_trunc <- ptrunc(
qt, "beta", shp1, shp2, a = a, lower.tail = lt, log.p = lg
)
p_beta <- pbeta(qt, shp1, shp2, ncp = 0, lt, lg)
if (!lg) {
expect_gte(p_trunc, 0)
expect_lte(p_trunc, 1)
if (lt) {
expect_lte(p_trunc, p_beta)
expect_length(qt, i)
expect_length(p_trunc, i)
for (q in seq_along(qt)) {
if (!lg) {
expect_gte(p_trunc[q], 0)
expect_lte(p_trunc[q], 1)
if (lt) {
expect_lte(p_trunc[q], p_beta[q])
} else {
expect_gt(p_trunc[q], p_beta[q])
}
} else {
expect_gt(p_trunc, p_beta)
expect_lte(p_trunc[q], 0)
}
} else {
expect_lte(p_trunc, 0)
}
}
}
Expand All @@ -63,20 +69,22 @@ test_that("lower truncation works as expected (beta)", {
test_that("lower truncation works as expected (binomial)", {
for (lt in c(TRUE, FALSE)) {
for (lg in c(FALSE, TRUE)) {
for (i in seq_len(10)) {
size <- sample(10:50, 1L)
prob <- runif(1)
for (i in seq_len(5)) {
size <- sample(10:30, 1L)
prob <- runif(1L, .2, .8)
a <- sample(1:(size - 4L), 1L)
qt <- sample(seq(a + 1L, size - 1L), 1L)
p_trunc <- ptrunc(
qt, "binomial", size, prob, a = a, lower.tail = lt, log.p = lg
)
qt <- sample(seq(a + 1L, size - 1L), i, replace = TRUE)
p_trunc <- ptrunc(qt, "binomial", size, prob, a = a, lower.tail = lt, log.p = lg)
p_binom <- pbinom(qt, size, prob, lower.tail = lt, log.p = lg)
if (!lg) {
expect_gte(p_trunc, 0)
expect_lte(p_trunc, 1)
} else {
expect_lte(p_trunc, 0)
expect_length(qt, i)
expect_length(p_trunc, i)
for (q in seq_along(qt)) {
if (!lg) {
expect_gte(p_trunc[q], 0)
expect_lte(p_trunc[q], 1)
} else {
expect_lte(p_trunc[q], 0)
}
}
}
}
Expand All @@ -86,20 +94,24 @@ test_that("lower truncation works as expected (binomial)", {
test_that("lower truncation works as expected (poisson)", {
for (lt in c(TRUE, FALSE)) {
for (lg in c(FALSE, TRUE)) {
for (i in seq_len(10)) {
for (i in seq_len(5)) {
lambda <- sample(10:50, 1L)
max_qt <- qpois(p = .99, lambda)
a <- sample(seq(1L, max_qt - 3L), 1L)
qt <- sample(seq(a + 1L, max_qt), 1L)
qt <- sample(seq(a + 1L, max_qt), i, replace = TRUE)
p_trunc <- ptrunc(
qt, "poisson", lambda, a = a, lower.tail = lt, log.p = lg
)
p_pois <- ppois(qt, lambda, lower.tail = lt, log.p = lg)
if (!lg) {
expect_gte(p_trunc, 0)
expect_lte(p_trunc, 1)
} else {
expect_lte(p_trunc, 0)
expect_length(qt, i)
expect_length(p_trunc, i)
for (q in seq_along(qt)) {
if (!lg) {
expect_gte(p_trunc[q], 0)
expect_lte(p_trunc[q], 1)
} else {
expect_lte(p_trunc[q], 0)
}
}
}
}
Expand All @@ -109,54 +121,155 @@ 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)) {
for (i in seq_len(5)) {
df <- sample(1:100, 1L)
a <- min(rchisq(10L, df))
qt <- max(rchisq(10L, df), a)
qt <- replicate(i, 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)
expect_length(qt, i)
expect_length(p_trunc, i)
for (q in seq_along(qt)) {
if (!lg) {
expect_gte(p_trunc[q], 0)
expect_lte(p_trunc[q], 1)
} else {
expect_lte(p_trunc[q], 0)
}
}
}
}
}
})

test_that("lower truncation works as expected (contbern)", {
for (i in seq_len(10)) {
for (i in seq_len(5)) {
lambda <- runif(1L)
a <- runif(1L)
qt <- runif(1L, a, 1L)
qt <- runif(i, a, 1L)
p_trunc <- ptrunc(qt, "contbern", lambda, a = a)
p_contbern <- pcontbern(qt, lambda)
expect_lte(p_trunc, p_contbern)
expect_length(qt, i)
expect_length(p_trunc, i)
for (q in seq_along(qt)) {
expect_lte(p_trunc[q], p_contbern[q])
}
}
})

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)) {
for (i in seq_len(5)) {
rate <- rchisq(1L, df = 10L)
a <- rexp(1L, rate)
qt <- max(rexp(10L, rate), a)
qt <- replicate(i, 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)
expect_length(qt, i)
expect_length(p_trunc, i)
for (q in seq_along(qt)) {
if (!lg) {
expect_gte(p_trunc[q], 0)
expect_lte(p_trunc[q], 1)
} else {
expect_lte(p_trunc[q], 0)

}
}
}
}
}
})

test_that("lower truncation works as expected (gamma)", {
for (lt in c(TRUE, FALSE)) {
for (lg in c(FALSE, TRUE)) {
for (i in seq_len(5)) {
shape <- rchisq(1L, df = 10L)
rate <- rchisq(1L, df = 10L)
a <- rgamma(1L, shape, rate)
qt <- replicate(i, max(rgamma(10L, shape, rate), a))
p_trunc <- ptrunc(
qt, "gamma", shape, rate, a = a, lower.tail = lt, log.p = lg
)
p_trunc_2 <- ptrunc(
qt, "gamma", shape, scale = 1 / rate, a = a, lower.tail = lt,
log.p = lg
)
p_gamma <- pgamma(qt, shape, rate, lower.tail = lt, log.p = lg)
expect_length(qt, i)
expect_length(p_trunc, i)
for (q in seq_along(qt)) {
if (!lg) {
expect_gte(p_trunc[q], 0)
expect_lte(p_trunc[q], 1)
expect_gte(p_trunc_2[q], 0)
expect_lte(p_trunc_2[q], 1)
} else {
expect_lte(p_trunc[q], 0)
expect_lte(p_trunc_2[q], 0)
}
}
expect_equal(p_trunc, p_trunc_2)
}
}
}
})

test_that("lower truncation works as expected (invgamma)", {
for (lt in c(TRUE, FALSE)) {
for (lg in c(FALSE, TRUE)) {
for (i in seq_len(5)) {
shape <- rchisq(1L, df = 10L)
rate <- rchisq(1L, df = 10L)
a <- rinvgamma(1L, shape, rate)
qt <- replicate(i, max(rinvgamma(10L, shape, rate), a))
p_trunc <- ptrunc(
qt, "invgamma", shape, rate, a = a, lower.tail = lt, log.p = lg
)
p_trunc_2 <- ptrunc(
qt, "invgamma", shape, scale = 1 / rate, a = a, lower.tail = lt,
log.p = lg
)
p_invgamma <- pinvgamma(qt, shape, rate, lower.tail = lt, log.p = lg)
expect_length(qt, i)
expect_length(p_trunc, i)
for (q in seq_along(qt)) {
if (!lg) {
expect_gte(p_trunc[q], 0)
expect_lte(p_trunc[q], 1)
expect_gte(p_trunc_2[q], 0)
expect_lte(p_trunc_2[q], 1)
} else {
expect_lte(p_trunc[q], 0)
expect_lte(p_trunc_2[q], 0)
}
}
expect_equal(p_trunc, p_trunc_2)
}
}
}
})

test_that("lower truncation works as expected (invgauss)", {
for (i in seq_len(5)) {
m <- rchisq(1L, df = 10L)
s <- rchisq(1L, df = 10L)
a <- rinvgauss(1L, m, s)
qt <- replicate(i, max(rinvgauss(10L, m, s), a))
p_trunc <- ptrunc(qt, "invgauss", m, s, a)
p_invgauss <- pinvgauss(qt, m, s)
expect_length(qt, i)
expect_length(p_trunc, i)
for (q in seq_along(qt)) {
expect_gte(p_trunc[q], 0)
expect_lte(p_trunc[q], 1)
expect_lte(p_trunc[q], p_invgauss[q])
}
}
})
Loading

0 comments on commit 5bb4d85

Please sign in to comment.