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
