Nothing
### This is a translation of R_zeroin2 in ~/R/D/r-devel/R/src/appl/zeroin.c
### from C to R by John Nash,
### ---> file rootoned/R/zeroin.R of the new (2011-08-18) R-forge package rootoned
###
### Where John Nash calls it zeroin(), I call it unirootR()
##' Simple modification of uniroot() which should work with mpfr-numbers
##' MM: uniroot() is in ~/R/D/r-devel/R/src/library/stats/R/nlm.R
##'
unirootR <- function(f, interval, ...,
lower = min(interval), upper = max(interval),
f.lower = f(lower, ...), f.upper = f(upper, ...),
extendInt = c("no", "yes", "downX", "upX"),
trace = 0,
verbose = as.logical(trace), verbDigits = max(3, min(20, -log10(tol)/2)),
tol = .Machine$double.eps^0.25, maxiter = 1000L,
check.conv = FALSE,
## Rmpfr-only:
warn.no.convergence = !check.conv,
epsC = NULL)
{
if(!missing(interval) && length(interval) != 2L)
stop("'interval' must be a vector of length 2")
## For many "quick things", we will use as.numeric(.) but we do *NOT* assume that
## lower and upper are numeric!
.N <- as.numeric
if(lower >= upper) # (may be mpfr-numbers *outside* double.xmax)
stop("lower < upper is not fulfilled")
if(is.na(.N(f.lower))) stop("f.lower = f(lower) is NA")
if(is.na(.N(f.upper))) stop("f.upper = f(upper) is NA")
form <- function(x, digits = verbDigits) format(x, digits=digits, drop0trailing=TRUE)
formI <- function(x, di = getOption("digits")) format(x, digits=di, drop0trailing=TRUE)
Sig <- switch(match.arg(extendInt),
"yes" = NULL,
"downX"= -1,
"no" = 0,
"upX" = 1,
stop("invalid 'extendInt'; please report"))
## protect against later 0 * Inf |--> NaN and Inf * -Inf.
truncate <- function(x) { ## NA are already excluded; deal with +/- Inf
if(is.numeric(x))
pmax.int(pmin(x, .Machine$double.xmax), -.Machine$double.xmax)
else if(inherits(x, "mpfr") && is.infinite(x)) # use maximal/minimal mpfr-number instead:
sign(x) * mpfr(2, .getPrec(x))^((1 - 2^-52)*.mpfr_erange("Emax"))
else x
}
f.low. <- truncate(f.lower)
f.upp. <- truncate(f.upper)
doX <- ( is.null(Sig) && f.low. * f.upp. > 0 ||
is.numeric(Sig) && (Sig*f.low. > 0 || Sig*f.upp. < 0))
if(doX) { ## extend the interval = [lower, upper]
if(trace)
cat(sprintf("search {extendInt=\"%s\", Sig=%s} in [%s,%s]%s",
extendInt, formI(Sig), formI(lower), formI(upper),
if(trace >= 2)"\n" else " ... "))
Delta <- function(u) 0.01* pmax(1e-4, abs(u)) ## <-- FIXME? [= R's uniroot() for double]
it <- 0L
## Two cases:
if(is.null(Sig)) {
## case 1) 'Sig' unspecified --> extend (lower, upper) at the same time
delta <- Delta(c(lower,upper))
while(isTRUE(f.lower*f.upper > 0) &&
any(iF <- is.finite(c(lower,upper)))) {
if((it <- it + 1L) > maxiter)
stop(gettextf("no sign change found in %d iterations", it-1),
domain=NA)
if(iF[1]) {
ol <- lower; of <- f.lower
if(is.na(f.lower <- f(lower <- lower - delta[1], ...))) {
lower <- ol; f.lower <- of; delta[1] <- delta[1]/4
}
}
if(iF[2]) {
ol <- upper; of <- f.upper
if(is.na(f.upper <- f(upper <- upper + delta[2], ...))) {
upper <- ol; f.upper <- of; delta[2] <- delta[2]/4
}
}
if(trace >= 2)
cat(sprintf(" .. modified lower,upper: (%15g,%15g)\n",
.N(lower), .N(upper)))
delta <- 2 * delta
}
} else {
## case 2) 'Sig' specified --> typically change only *one* of lower, upper
## make sure we have Sig*f(lower) <= 0 and Sig*f(upper) >= 0:
delta <- Delta(lower)
while(isTRUE(Sig*f.lower > 0)) {
if((it <- it + 1L) > maxiter)
stop(gettextf("no sign change found in %d iterations", it-1),
domain=NA)
f.lower <- f(lower <- lower - delta, ...)
if(trace >= 2) cat(sprintf(" .. modified lower: %s, f(.)=%s\n",
formI(lower), formI(f.lower)))
delta <- 2 * delta
}
delta <- Delta(upper)
while(isTRUE(Sig*f.upper < 0)) {
if((it <- it + 1L) > maxiter)
stop(gettextf("no sign change found in %d iterations", it-1),
domain=NA)
f.upper <- f(upper <- upper + delta, ...)
if(trace >= 2) cat(sprintf(" .. modified upper: %s, f(.)=%s\n",
formI(upper), formI(f.upper)))
delta <- 2 * delta
}
}
if(trace && trace < 2)
cat(sprintf("extended to [%s, %s] in %d steps\n", formI(lower), formI(upper), it))
}
if(!isTRUE(sign(f.lower) * sign(f.upper) <= 0))
stop(if(doX)
"did not succeed extending the interval endpoints for f(lower) * f(upper) <= 0"
else sprintf("f() values at end points = (%s, %s) not of opposite sign",
formI(f.lower), formI(f.upper)))
if(is.null(epsC) || is.na(epsC)) {
## determine 'epsC' ``the achievable Machine precision''
## -- given the class of f.lower, f.upper
ff <- f.lower * f.upper
if(is.double(ff)) epsC <- .Machine$double.eps
else if(is(ff, "mpfr"))
epsC <- 2^-min(getPrec(f.lower), getPrec(f.upper))
else { ## another number class -- try to see if getPrec() is defined..
## if not, there's not much we can do
if(is(prec <- tryCatch(min(getPrec(f.lower), getPrec(f.upper)),
error = function(e)e),
"error")) {
warning("no valid getPrec() for the number class(es) ",
paste(unique(class(f.lower),class(f.upper)), collapse=", "),
".\n Using double precision .Machine$double.eps.")
epsC <- .Machine$double.eps
} else {
epsC <- 2^-prec
message("using epsC = %s ..", format(epsC))
}
}
}
if(tol < epsC / 8) # "8 fudge factor" (otherwise happens too often)
warning(sprintf("tol (%g) < epsC (%g) is rarely sensical,
and the resulting precision is probably not better than epsC",
tol, epsC))
## Instead of the call to C code, now "do it in R" :
## val <- .Internal(zeroin2(function(arg) as.numeric(f(arg, ...)),
## lower, upper, f.lower, f.upper,
## tol, as.integer(maxiter)))
a <- lower # interval[1]
b <- upper # interval[2]
fa <- f.lower # f(ax, ...)
fb <- f.upper # f(bx, ...)
if (verbose)
cat(sprintf("==> Start zeroin: f(%g)= %g; f(%g)= %g\n", .N(a), .N(fa), .N(b), .N(fb)))
c <- a
fc <- fa
## First test if we have found a root at an endpoint
maxit <- maxiter + 2L # count evaluations as maxiter-maxit
converged <- FALSE
while(!converged && maxit > 0) { ##---- Main iteration loop ------------------------------
if (verbose) cat("Iteration >>>", maxiter+3L-maxit, "<<< ;")
d.prev <- b-a
## Distance from the last but one to the last approximation */
##double tol.2; ## Actual tolerance */
##double p; ## Interpolation step is calcu- */
##double q; ## lated in the form p/q; divi-
## * sion operations is delayed
## * until the last moment */
##double d.new; ## Step at this iteration */
if(abs(fc) < abs(fb)) { ## Swap data for b to be the smaller
if (verbose) cat(sprintf("fc (=%s) smaller than fb\n", form(fa)))
a <- b
b <- c
c <- a ## best approximation
fa <- fb
fb <- fc
fc <- fa
}
tol.2 <- 2*epsC*abs(b) + tol/2
d.new <- (c-b)/2 # bisection
if (verbose) cat("tol.2(epsC,b) = ",.N(tol.2), "; d.new= ",.N(d.new),"\n", sep="")
## converged <- (abs(d.new) <= tol.2 && is.finite(fb)) || fb == 0
converged <- (abs(d.new) <= tol.2) || fb == 0
if(converged) {
if (verbose) cat("DONE! -- small d.new or fb=0\n")
## Acceptable approx. is found :
val <- list(root=b, froot=fb, rtol = abs(c-b), maxit=maxiter-maxit)
}
else {
## Decide if the interpolation can be tried */
if( (abs(d.prev) >= tol.2) ## If d.prev was large enough*/
&& (abs(fa) > abs(fb)) ) { ## and was in true direction,
## Interpolation may be tried */
## register double t1,cb,t2;
if (verbose) cat("d.prev larger than tol.2 and fa bigger than fb --> ")
cb <- c-b
if (a == c) { ## If we have only two distinct points, linear interpolation
## can only be applied
t1 <- fb/fa
p <- cb*t1
q <- 1 - t1
if (verbose) cat("a == c: ")
}
else { ## Quadric inverse interpolation*/
if (verbose) cat("a != c: ")
q <- fa/fc
t1 <- fb/fc
t2 <- fb/fa
p <- t2 * ( cb*q*(q-t1) - (b-a)*(t1-1) )
q <- (q-1) * (t1-1) * (t2-1)
}
if(p > 0) { ## p was calculated with the */
if (verbose) cat(" p > 0; ")
q <- -q ## opposite sign; make p positive */
} else { ## and assign possible minus to */
if (verbose) cat(" p <= 0; ")
p <- -p ## q */
}
if (p < 0.75*cb*q - abs(tol.2*q)/2 ## If b+p/q falls in [b,c]*/
&& p < abs(d.prev*q/2)) { ## and isn't too large */
if (verbose) cat("p satisfies conditions for changing d.new\n")
d.new <- p/q ## it is accepted
} else if(verbose) cat("\n")
## If p/q is too large, then the
## bisection procedure can reduce [b,c] range to more extent
}
if( abs(d.new) < tol.2) { ## Adjust the step to be not less than tolerance
if (verbose) cat("d.new smaller than tol.2, adjusted to it.\n")
d.new <- if(d.new > 0) tol.2 else -tol.2
}
a <- b
fa <- fb ## Save the previous approx. */
b <- b + d.new
fb <- f(b, ...)
if (verbose) cat(sprintf("new f(b=%s) = %s;\n", form(b), form(fb)))
maxit <- maxit-1
## Do step to a new approxim. */
if( ((fb > 0) && (fc > 0)) || ((fb < 0) && (fc < 0)) ) {
if (verbose) cat(sprintf(" make c:=a=%s to have sign opposite to b: f(c)=%s\n",
form(a), form(fa)))
## Adjust c for it to have a sign opposite to that of b */
c <- a
fc <- fa
}
}## else not converged
} ## end{ while(maxit > 0) } --------------------------------------------
if(converged) {
iter <- val[["maxit"]]
if(!is.na(fb) && abs(fb) > 0.5*max(abs(f.lower), abs(f.upper)))# from John Nash:
warning("Final function magnitude seems large -- maybe converged to sign-changing 'pole' location?")
} else { ## (!converged) : failed!
if(check.conv)
stop("no convergence in zero finding in ", iter, " iterations")
## else
val <- list(root= b, rtol = abs(c-b))
iter <- maxiter
if(warn.no.convergence)
warning("_NOT_ converged in ", iter, " iterations")
}
list(root = val[["root"]], f.root = f(val[["root"]], ...),
iter = iter, estim.prec = .N(val[["rtol"]]), converged = converged)
} ## {unirootR}
### qnorm() via inversion of pnorm() by unirootR() ==========================
"
qnorm(p) == q <==> pnorm(q) == p
"
##====== TODO: Use good i.e. *tight* inequalities for Phi(x) .. which are invertible
## ===== to (still tight) inequalities for Phi^{-1}(x) == qnorm(x)
## ===> get good *start* interval for unirootR()
## qnorm (p, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
qnormI <- function(p, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE,
trace = 0, verbose = as.logical(trace), # <- identical defaults as unirootR()
tol, # for "base" = .Machine$double.eps^0.25, but here use getPrec()
useMpfr = any(prec > 53), # if true, use mpfr
give.full = FALSE,
...) # <- arguments to pass to unirootR()
{
## The function whose "zeros" aka "roots" we want to find:
zFun <- function(q) pnorm(q, mean=mean, sd=sd, lower.tail=lower.tail, log.p=log.p) - p.
if(missing(tol) || !is.finite(tol)) {
prec <- max(getPrec(if(missing(mean) && missing(sd)) p else p+(mean+sd)))
## not max(getPrec(c(p, mean, sd))) as it's >= 128 by default
tol <- 2^-(prec + 2) # 2^(-prec - 1) gives less accurate
} else
prec <- as.integer(ceiling(1 - log2(tol)))
##
if(verbose) cat(sprintf("prec=%d ==> useMpfr=%s:\n", prec, format(useMpfr)))
verbDigs <- if(any(is.vd <- "verbDigits" == ...names()))
...elt(which(is.vd))
else max(3, min(20, -log10(tol)/2))
if(useMpfr) { ## This **IS** important here:
old_eranges <- .mpfr_erange() # typically -/+ 2^30
myERng <- (1-2^-52) * .mpfr_erange(c("min.emin","max.emax"))
if(!isTRUE(all.equal(myERng, old_eranges))) {
.mpfr_erange_set(value = myERng)
on.exit( .mpfr_erange_set( , old_eranges) )
}
}
sgn <- if(lower.tail) -1 else 1
INf <- if(lower.tail) Inf else -Inf
## "start"-interval for unirootR() --- see 'TODO' .. *tight* .. above <<<<<<<<
## Notably for the relevant (p <<< -1, log.p=TRUE) case !
qnInt <-
if(log.p) {
Pi <- if(useMpfr) Const("pi", prec) else pi
function(p) { s2 <- -2*p
if(useMpfr && !inherits(s2, "mpfr")) s2 <- mpfr(s2, precBits = prec)
xs1 <- s2 - log(2*Pi*s2)
if(p < -54) {
qn <- sqrt(s2 - log(2*Pi*xs1)); e <- 1e-4
} else if(p <= -15) {
qn <- sqrt(s2 - log(2*Pi*xs1) - 1/(2 + xs1)); e <- 1e-3
} else { # p >= -15
qn <- stats__qnorm(asNumeric(p), log.p=TRUE)
## FIXME: not good enough, e.g., for p = mpfr(-1e-5, 128)
e <- 1e-2
}
qn*(sgn + c(-e,e))
}
} else { ## log.p is FALSE
Id <- if(useMpfr) function(.) mpfr(., precBits = prec) else identity
function(p) { q <- stats__qnorm(asNumeric(p), lower.tail=lower.tail)
Id(if(abs(q) < 1e-3) c(-2,2)*1e-3 else c(.99, 1.01) * q)
}
}
## Deal with prob in {0, 1} which correspond to quantiles -/+ Inf :
## idea from {DPQ}'s .D_0 and .D_1 :
.p_1 <- as.integer(!log.p)
.p_0 <- if (log.p) -Inf else 0
r <- if(give.full) vector("list", length(p))
else if(useMpfr) mpfr(p, precBits=prec)
else p
for(ip in seq_along(p)) {
p. <- p[ip]
if(verbose) cat(sprintf("p. = p[ip=%d] = %s:\n", ip,
format(p., digits = verbDigs, drop0trailing=TRUE)))
ri <-
if (p. == .p_1)
INf
else if (p. == .p_0)
-INf
else if(is.na(p.) || p. > .p_1 || p. < .p_0)
NaN
else { ## zFun() uses p.
ur <- unirootR(zFun, interval = qnInt(p.),
extendInt = if(lower.tail) "upX" else "downX",
trace=trace, verbose=verbose, tol=tol,
...) # verbDigits, maxiter, check.conv, warn.no.convergence, epsC
if(give.full) ur else ur$root
}
if(give.full) r[[ip]] <- ri else r[ip] <- ri
}
r
}
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.