Skip to content

Commit

Permalink
Fixed buggy behavior when alpha = 0
Browse files Browse the repository at this point in the history
  • Loading branch information
marberts committed Feb 22, 2024
1 parent 0b3e743 commit a4569eb
Show file tree
Hide file tree
Showing 4 changed files with 21 additions and 11 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Package: sps
Title: Sequential Poisson Sampling
Version: 0.5.3.9008
Version: 0.5.3.9009
Authors@R: c(
person("Steve", "Martin",
role = c("aut", "cre", "cph"),
Expand Down
4 changes: 4 additions & 0 deletions R/inclusion_prob.R
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,9 @@ pi <- function(x, n, alpha, cutoff) {
stop("'cutoff' must be greater than 0")
}

# Add some fuzz to avoid floating-point rounding stopping some units from
# being TA, especially when alpha = 0.
alpha <- alpha + .Machine$double.eps^0.5
ta <- which(x >= cutoff)
if (length(ta) > n) {
stop("'n' is not large enough to include all units with 'x' above 'cutoff'")
Expand Down Expand Up @@ -181,6 +184,7 @@ becomes_ta <- function(x, alpha = 1e-3, cutoff = Inf) {
stop("'cutoff' must be greater than 0")
}

alpha <- alpha + .Machine$double.eps^0.5
ta <- which(x >= cutoff)
x[ta] <- 0
ord <- rev(order(x, decreasing = TRUE))
Expand Down
2 changes: 1 addition & 1 deletion R/sps.R
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ order_sampling_ <- function(f) {
ta <- which(ta)
n_ts <- n - length(ta)
# sample the take somes
keep <- if (n_ts > 0) {
keep <- if (n_ts > 0L) {
xi <- f(u[ts]) / f(p[ts])
order(xi)[seq_len(n_ts)]
}
Expand Down
24 changes: 15 additions & 9 deletions tests/testthat/test-inclusion_prob.R
Original file line number Diff line number Diff line change
@@ -1,15 +1,12 @@
set.seed(14235)

test_that("corner cases work as expected", {
expect_equal(
inclusion_prob(0, 0),
0
)
expect_equal(inclusion_prob(0, 0), 0)
expect_equal(
inclusion_prob(1:3, c(0, 1, 0), factor(c(2, 2, 2), levels = 1:3)),
1:3 / 6
)
expect_equal(
expect_identical(
inclusion_prob(1:6, c(0, 3), c(1, 1, 2, 1, 2, 2)),
c(0, 0, 1, 0, 1, 1)
)
Expand All @@ -21,14 +18,18 @@ test_that("corner cases work as expected", {
inclusion_prob(rep(1, 6), c(2, 1), c(1, 1, 2, 1, 2, 2)),
c(2, 2, 1, 2, 1, 1) / 3
)
expect_equal(
expect_identical(
inclusion_prob(c(0, 1, 1, 1 + 1e-4), 3),
c(0, 1, 1, 1)
)
expect_equal(
expect_identical(
inclusion_prob(c(0, 1, 1, 1 - 1e-4), 3),
c(0, 1, 1, 1)
)
expect_identical(
inclusion_prob(rep(0.1, 3), 3, alpha = 0),
c(1, 1, 1)
)

expect_equal(becomes_ta(0), NaN)
expect_equal(becomes_ta(1), 1)
Expand Down Expand Up @@ -212,6 +213,11 @@ test_that("adding a cutoff just offsets when a unit becomes TA", {
expect_equal(becomes_ta(x, 0.25), c(3, 4, 5, 4, 5))
expect_equal(becomes_ta(x, 0.25, 4), c(NaN, 4, 5, NaN, 5))

x <- c(10, 10, 10, x)
expect_equal(becomes_ta(x, 0.25, 6), c(NaN, NaN, NaN, 6, 7, 8, 7, 8))
y <- c(10, 10, 10, x)
expect_equal(becomes_ta(y, 0.25, 6), c(NaN, NaN, NaN, 6, 7, 8, 7, 8))

# Adding zeroes does nothing
z <- c(0, 0, 0, x)
expect_equal(becomes_ta(z, 0.25, 6),
c(NaN, NaN, NaN, 3, 4, 5, 4, 5))
})

0 comments on commit a4569eb

Please sign in to comment.