Nothing
### Was part of ./qnbinom-slow-ex.R
### ~~~~~~~~~~~~~~~~~
## MM_FIXME_: solve "small 'size'" cases: do_search() is really "too cheap" sometimes, there
## ---
if(FALSE) ## C source:
invisible( "~/R/D/r-devel/R/src/nmath/qbinom.c" )
##======= TODO: this is practically identical to doSearch_pnb() in ./qnbinomR.R === use ONE function
## and doSearch_poi() in ./qpoisR.R ... " " "
## with (<pDist>, <pDname>) arguments, here (pnbinom, "pnbinom")
## ==> can MPFR-ify it by passing pnbinomMpfr()..
##' @title For P(x) := pbinom(x, n, pr), find p-quantile y = y(p) :<==> P(y-1) < p <= P(y)
##' @param y current guess
##' @param z == pbinom(y, ..)
##' @param p target probability
##' @param n, pr : (n = size); parameters of the binomial
##' @param incr increment, integer valued >= 1.
##' @param lower.tail
##' @param log.p
##' @param trace
##' @return root 'y' and z[0] = pbinom(y, ..)
##' @author Martin Maechler
doSearch_bin <- function(y, z, p, n, pr, incr,
lower.tail=TRUE, log.p=FALSE, trace = 0)
{
formatI <- function(x) format(x, scientific = 16)
if(y < 0) y <- 0
iStr <- if(incr != 1) sprintf(" incr = %s", formatI(incr)) else ""
left <- if(lower.tail) (z >= p) else (z < p)
it <- 0L
if(left) { ## search to the left
if(trace) cat(sprintf("doSearch() to the left (pbinom(y=%g, *) = z = %g %s p=%g):%s\n",
y,z, if(lower.tail) ">=" else "<", p, iStr))
repeat {
it <- it + 1L
if(y > 0)
newz <- pbinom(y - incr, n, pr, lower.tail=lower.tail, log.p=log.p)
else if(y < 0)
y <- 0
if(y == 0 || is.na(newz) || ## pbinom() gives NaN rarely because of pbeta(*, log.p=TRUE) bug
if(lower.tail) (newz < p) else (newz >= p)) {
if(trace >= 2) cat(sprintf(
" new y=%.15g, pbinom(y-incr,*) %s p; ==> doSearch() returns prev z=%g after %d iter.\n",
y, if(lower.tail) "<" else ">=", z, it))
## pbinom(y-incr, *) < p <= prev.z = pbinom(y, *)
return(list(y=y, bin.y = z, iter = it))
}
y <- max(0, y - incr)
z <- newz
}
}
else { ## !left : search to the right
if(trace) cat(sprintf("doSearch() to the right (pbinom(y=%g, *) = z = %g %s p=%g):%s\n",
y,z, if(lower.tail) "<" else ">=", p, iStr))
repeat {
it <- it + 1L
y <- y + incr
z <- pbinom(y, n, pr, lower.tail=lower.tail, log.p=log.p)
if(is.na(z) || ## pbinom() returning NaN {pbeta() bug}
if(lower.tail) z >= p else z < p) {
if(trace >= 2) cat(sprintf(
" new y=%.15g, z=%g = pbinom(y,*) %s p ==> doSearch() returns after %d iter.\n",
y, z, if(lower.tail) ">=" else "<", it))
## pbinom(y, *) >= p > prev.z = pbinom(prev.y, *)= pbinom(y-incr, *) <===>
## pbinom(y-incr, *) = prev.z < p <= pbinom(y, *)
return(list(y=y, bin.y = z, iter = it))
}
}
}
}
qbinomR1 <- function(p, size, prob, lower.tail=TRUE, log.p=FALSE,
yLarge = 4096, # was hard wired to 1e5
incF = 1/64, # was hard wired to .001
iShrink = 8, # was hard wired to 100
relTol = 1e-15,# was hard wired to 1e-15
pfEps.n = 8, # was hard wired to 64: "fuzz to ensure left continuity"
pfEps.L = 2, # was hard wired to 64: " " ..
fpf = 4, # *MUST* be >= 1 (did not exist previously)
trace = 0)
{
stopifnot(length(p) == 1, length(size) == 1, length(prob) == 1)
# #ifdef IEEE_754
if (is.na(p) || is.na(size) || is.na(prob))
return(p + size + prob)
# #endif
if (prob <= 0 || prob > 1 || size < 0) {
warning("returning NaN because of (prob,size) out of range")
return(NaN)
}
if (prob == 1 || size == 0) return(0)
## R_Q_P01_boundaries(p, 0, size); (n == size)
## {TODO? simplify via if(p == .D_0(log.p)) .. as in ./t-dist.R }
if(log.p) {
if(p == -Inf) return(if(lower.tail) 0 else size)
if(p == 0 ) return(if(lower.tail) size else 0)
if(p > 0) { warning("p > 0 -> returning NaN"); return(NaN) }
## now p \in (-Inf, 0) (open interval)
} else {
if(p == 0) return(if(lower.tail) 0 else size)
if(p == 1) return(if(lower.tail) size else 0)
if(p < 0 || p > 1) {
warning("p out of [0,1] -> returning NaN"); return(NaN) }
## now p \in (0, 1) (open interval)
}
q <- 1 - prob;
mu <- size * prob;
sigma <- sqrt(size * prob * q);
gamma <- (q - prob) / sigma;
if(trace) {
cat(sprintf("qbinom(p=%.12g, size=%.15g, prob=%g, l.t.=%d, log=%d):\n
mu=%g, sigma=%g, gamma=%g;\n",
p, size, prob, lower.tail, log.p, mu, sigma, gamma))
}
## /* Note : "same" code in qpois.c, qbinom.c, qbinom.c --
## * FIXME: This is far from optimal [cancellation for p ~= 1, etc]: */
if(!lower.tail || log.p) {
p_n <- .DT_qIv(p, lower.tail, log.p) ## need check again (cancellation!):
if(trace) cat(sprintf(" upper tail or log_p: (p_n, 1-p_n) = (%.15g, %.15g)\n", p_n, 1-p_n));
if (p_n == 0) message("p_n=0: NO LONGER return(0)\n")
if (p_n == 1) message("p_n=1: NO LONGER return(Inf)\n")
} else
p_n <- p
## temporary hack --- FIXME --- */
if (p_n + 1.01*.Machine$double.eps >= 1.) {
if(trace)
cat("p_n + 1.01 * c_eps >= 1 ; (1-p = ",format(1-p),
") ==> NOW LONGER returning Inf (Hack -- FIXME ?)\n", sep="")
## return(Inf)
}
## y := approx.value (Cornish-Fisher expansion) : */
z <- qnorm(p, lower.tail=lower.tail, log.p=log.p)
y <- round(mu + sigma * (z + gamma * (z*z - 1) / 6)) # R_forceint(.)
if(trace) cat(sprintf("Cornish-Fisher: initial z=%g, y=%g; ", z,y))
if(y > size) # way off
y <- size
else if(y < 0) ## (new check; (when) does it happen ?? )
y <- 0.
z <- pbinom(y, size, prob, lower.tail=lower.tail, log.p=log.p)
if(trace) cat(sprintf(" then: y=%g, z= pbinom(y,*) = %g;", y,z))
## fuzz to ensure left continuity: */
##p <- p*( 1 - pEps.f *.Machine$double.eps ) ## pEps.f was hard wired to 64
## fuzz to ensure left continuity: do not loose too much (=> error in upper tail)
c.eps <- .Machine$double.eps
if(log.p) { ## <==> p \in [-Inf, 0] different adjustment: "other sign"
e <- pfEps.L * c.eps
if(lower.tail && p > -.Machine$x.max)# prevent underflow to -Inf
p <- p * (1 + e)
else ## if(p < - .Machine$x.min) # not too close to -0
p <- p * (1 - e)
} else { ##// not log scale
e <- pfEps.n * c.eps
if(lower.tail)
p <- p * (1 - e)
else if(1 - p > fpf*e)# otherwise get p > 1
p <- p * (1 + e)
}
if(trace) cat(sprintf(" left-cont. fuzz => p = %15.11g\n", p))
## If the C-F value is not too large a simple search is OK */
if(y < yLarge) {
doSearch_bin(y=y, z=z, p=p, n=size, pr=prob, incr = 1,
lower.tail=lower.tail, log.p=log.p, trace=trace)$ y
## pbinom(y-1, *) < p <= pbinom(y, *) ==> y is solution
} else { ## y >= yLarge, be a bit cleverer in the search:
## use larger increments, notably initially:
incr <- floor(y * incF) ## incF was hard-wired to 0.001
if(trace) cat(sprintf("large y: --> use larger increments than 1: incr=%s\n",
format(incr, scientific=16)))
qIt <- totIt <- 0L
repeat {
oldincr <- incr
yz <- doSearch_bin(y=y, z=z, p=p, n=size, pr=prob, incr = incr,
lower.tail=lower.tail, log.p=log.p, trace=trace)
## We know (for both, left or right search !) :
## pbinom(y-oldincr, *) < p <= pbinom(y, *)
y <- yz$y
z <- yz$bin.y
totIt <- totIt + yz$iter
qIt <- qIt + 1L
incr <- max(1, floor(incr/iShrink)); ## was (incr/100)
## while(oldincr > 1 && incr > y*1e-15);
if(oldincr <= 1 || incr <= y * relTol)
break
}
if(trace) cat(sprintf(" \\--> %s; needed %d doSearch() calls, total %d iterations\n",
if(oldincr == 1) "oldincr=1"
else sprintf("oldincr = %s, incr = %s < %.11g = y * relTol",
format(oldincr, scientific = 9),
format(incr, scientific = 9), y*relTol),
qIt, totIt))
## return
y
}
} ## qbinomR1()
## A version vectorized in 'p' :
qbinomR <- Vectorize(qbinomR1, c("p", "size", "prob"))
# but not {lower.tail, log.p, iShrink, ....}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.