>>>>> Ravi Varadhan via R-devel writes:
Thanks. I finally committed
r90155 | hornik | 2026-06-15 18:45:43 +0200 (Mon, 15 Jun 2026) | 2 lines
Add experimental code for faster computation of minlike two-sided
p-values for large n. Based on R-devel post by Ravi Varadhan.
Needed to make a few other changes first. I want to experiment some
more in the next few days and then fold in, likely already for n around
1000 (for which improvements become noticeable for me).
One should also modify the code in poisson.test() accordingly ...
-k
> This code block in binom.test() is the bottleneck for large n:
> else if (x < m) {
> i <- seq.int(from = ceiling(m), to = n)
> y <- sum(dbinom(i, n, p) <= d * relErr)
> pbinom(x, n, p) +
> pbinom(n - y, n, p, lower.tail = FALSE)
> } else {
> i <- seq.int(from = 0, to = floor(m))
> y <- sum(dbinom(i, n, p) <= d * relErr)
> pbinom(y - 1, n, p) +
> pbinom(x - 1, n, p, lower.tail = FALSE)
> }
> Instead of summing over the entire support, we can:
> 1.
> Find the left and right boundary integers, yL and yR, where the PMF drops to
> or below the observed probability (using uniroot).
> 2.
> Compute the exact p-value as: pbinom(yL, n, p) + (1 - pbinom(yR-1, n, p)).
> This is much faster: O(log n) for the search versus O(n).
> Here is a function (I used n=1e05 as the threshold):
> binom.test.fast <- function(x, n, p = 0.5, alternative = c("two.sided",
> "less", "greater"), conf.level = 0.95) {
> DNAME <- deparse1(substitute(x))
> xr <- round(x)
> if (any(is.na(x) | (x < 0)) || max(abs(x - xr)) > 1e-07)
> stop("'x' must be nonnegative and integer")
> x <- xr
> if (length(x) == 2L) {
> n <- sum(x)
> x <- x[1L]
> }
> else if (length(x) == 1L) {
> nr <- round(n)
> if ((length(n) > 1L) || is.na(n) || (n < 1) || abs(n - nr) > 1e-07 || (x
> > nr))
> stop("'n' must be a positive integer >= 'x'")
> DNAME <- paste(DNAME, "and", deparse1(substitute(n)))
> n <- nr
> }
> else stop("incorrect length of 'x'")
> if (!missing(p) && (length(p) > 1L || is.na(p) || p < 0 || p > 1))
> stop("'p' must be a single number between 0 and 1")
> alternative <- match.arg(alternative)
> if (!((length(conf.level) == 1L) && is.finite(conf.level) && (conf.level >
> 0) && (conf.level < 1)))
> stop("'conf.level' must be a single number between 0 and 1")
> PVAL <- switch(alternative,
> less = pbinom(x, n, p),
> greater = pbinom(x - 1, n, p, lower.tail = FALSE),
> two.sided = {
> if (p == 0) (x == 0) else if (p == 1) (x == n) else {
> if (n > 1e5) {
> # Fast O(log n) path for large n
> log_dx <- dbinom(x, n, p, log = TRUE)
> if (!is.finite(log_dx)) {
> 0
> } else {
> relErr <- 1 + 1e-07
> log_thresh <- log_dx + log(relErr)
> mode <- floor((n + 1) * p)
> # Left boundary: largest k <= mode with dbinom(k) <=
> dbinom(x)*relErr
> lo <- 0L; hi <- min(x, mode)
> while (lo < hi) {
> mid <- (lo + hi + 1L) %/% 2L
> if (dbinom(mid, n, p, log = TRUE) > log_thresh) hi
> <- mid - 1L
> else lo <- mid
> }
> yL <- lo
> # Right boundary: smallest k >= mode with dbinom(k)
> <= dbinom(x)*relErr
> lo <- max(x, mode); hi <- n
> while (lo < hi) {
> mid <- (lo + hi) %/% 2L
> if (dbinom(mid, n, p, log = TRUE) > log_thresh) lo
> <- mid + 1L
> else hi <- mid
> }
> yR <- lo
> pbinom(yL, n, p) + pbinom(yR, n, p, lower.tail =
> FALSE)
> }
> } else {
> # Original vectorized path for small n
> relErr <- 1 + 1e-07
> d <- dbinom(x, n, p)
> m <- n * p
> if (x == m) 1 else if (x < m) {
> i <- seq.int(from = ceiling(m), to = n)
> y <- sum(dbinom(i, n, p) <= d * relErr)
> pbinom(x, n, p) + pbinom(n - y, n, p, lower.tail =
> FALSE)
> } else {
> i <- seq.int(from = 0, to = floor(m))
> y <- sum(dbinom(i, n, p) <= d * relErr)
> pbinom(y - 1, n, p) + pbinom(x - 1, n, p, lower.tail
> = FALSE)
> }
> }
> }
> }
> )
> p.L <- function(x, alpha) {
> if (x == 0) 0 else qbeta(alpha, x, n - x + 1)
> }
> p.U <- function(x, alpha) {
> if (x == n) 1 else qbeta(1 - alpha, x + 1, n - x)
> }
> CINT <- switch(alternative,
> less = c(0, p.U(x, 1 - conf.level)),
> greater = c(p.L(x, 1 - conf.level), 1),
> two.sided = {
> alpha <- (1 - conf.level)/2
> c(p.L(x, alpha), p.U(x, alpha))
> }
> )
> attr(CINT, "conf.level") <- conf.level
> ESTIMATE <- x/n
> names(x) <- "number of successes"
> names(n) <- "number of trials"
> names(ESTIMATE) <- names(p) <- "probability of success"
> structure(list(statistic = x, parameter = n, p.value = PVAL,
> conf.int = CINT, estimate = ESTIMATE, null.value = p,
> alternative = alternative, method = "Exact binomial test",
> data.name = DNAME), class = "htest")
> }
> Ravi
> ________________________________
> From: R-devel <[email protected]> on behalf of Martin Maechler
> <[email protected]>
> Sent: Monday, May 18, 2026 12:20
> To: Chris Chang <[email protected]>
> Cc: R Development List <[email protected]>
> Subject: Re: [Rd] Feature request: ..., binomial ... exact test efficiency
> improvements
> External Email - Use Caution
>>>>> Chris Chang
>>>>> on Sat, 16 May 2026 14:24:19 -0700 writes:
> [............................]
>>> The one-sided binom.test() does not have this kind of
>>> performance problem, but the two-sided test (which can be
>>> written in a manner similar to qbinom(): instead of
>>> searching for minimal k satisfying pbinom(k, ...) >= p,
>>> search for the innermost opposite-tail k satisfying
>>> dbinom(k) <= dbinom(x)) does.
>>>
>>> I can submit a patch with efficient two-sided binomial
>>> and Fisher's 2x2 exact test implementations if there is
>>> interest.
>>>
>>> --Chris
> Dear Chris (and other readers),
> Above, you also mention binom.test(),
> and indeed, when I choose numbers in the order of 1e8,
> I see binom.test() taking quite some time, e.g.,
>> system.time(btL <- binom.test(c(4e8, 1.01e8), p = 0.8, alternative =
>> "two.sided"))
> user system elapsed
> 8.306 0.837 9.159
> and here, the time spent is really during p-value computation itself,
> everything else in binom.test() being very fast, AFAICS.
> Here it would be useful if you provide a patch to speed the
> p-value computation, still getting exact p-values, using
> pbinom() only instead of dbinom(<large_support), ..).
> Traditionally, we as R core and any applied statistician would
> say that in such large N cases, simple normal approximations should be
> accurate enough,
> which it would be practically in any case, OTOH, it would still
> be nice to get the exact probabilities / confidence intervals,
> with a faster calulation.
> Eventually using a new argument `exact = N < 1e6` (for some N)
> and using a normal (or better!) asymptotic approximation for
> exact = FALSE may make sense here as well...
> Thank you for raising the issue!
> With best regards,
> Martin
> --
> Martin Maechler
> ETH Zurich and R Core team
> ______________________________________________
> [email protected] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel<https://stat.ethz.ch/mailman/listinfo/r-devel>
> [[alternative HTML version deleted]]
> ______________________________________________
> [email protected] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
