R devel - May 2026 - Feature request: ..., binomial ... exact test efficiency improvements

>>>>> 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

>>>>> Martin Maechler writes:
Hmm.  binom.test is documented to perform an "Exact Binomial Test"
...?

-k
>>>>> 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
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

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 <r-devel-bounces at r-project.org> on behalf of Martin
Maechler <maechler at stat.math.ethz.ch>
Sent: Monday, May 18, 2026 12:20
To: Chris Chang <chrchang523 at gmail.com>
Cc: R Development List <R-devel at r-project.org>
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

______________________________________________
R-devel at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel<https://stat.ethz.ch/mailman/listinfo/r-devel>

	[[alternative HTML version deleted]]

Hi,

In the binomial and Fisher?s 2x2 cases, the exact p-value calculation
remains tractable for N well beyond 1e6.  Ravi's function shows roughly
what I had in mind for binom.test, and the fisher.test two-sided p-value
calculation (as well as other two-sided exact tests based on unimodal
distributions with an efficient cdf+ccdf implementation) would work
essentially the same way.

(In the meantime, I'm drafting https://pypi.org/project/exact-tests/ ,
which uses a "double-double" arithmetic library to mitigate the
log-factorial catastrophic cancellation problem.  It's still a few weeks
away from what I'd consider production quality, but it does already contain
accurate and efficient implementations of binomial and Fisher 2x2 exact
tests allowing N up to ~2^31 in general, and ~2^52 in the one-sided cases.
That is where I will initially try to implement efficient alternatives to
the current fisher.test odds-ratio and confidence-interval calculations.)

--Chris

On Mon, May 18, 2026 at 9:20?AM Martin Maechler <maechler at
stat.math.ethz.ch>
wrote:
> >>>>> 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
>
>
>
>
	[[alternative HTML version deleted]]