Nothing
require(DPQmpfr)
require(Rmpfr)
## very accurate set of *large* abscissa values
x <- mpfr(2,1024)^ seq(3, 100, by=1/2)
xN <- asNumeric(x)
## Looking at pnormU_S53() = upper bound to pnorm(), for large x,
## the term log(4 / (3 + sqrt(1 + 8/x^2))) early underflows to 0:
eps_S53 <- function(x, method = c("direct", "log1pExct", "deltaFull",
"delta1", "delta2", "delta3")) {
method <- match.arg(method)
delta <- function(x) { xm2 <- x^-2; (.5*(1 - sqrt(1 + 8*xm2)) + xm2) / (1 + xm2) }
switch(method,
"direct" = log(4 / (3 + sqrt(1 + 8*x^-2))),
"log1pExct" = log1p(delta(x)),
"deltaFull" = delta(x),
"delta1" = -1/(x^2+1) , ## 1 term for delta(x)
"delta2" = -(1 - (2/x)^2) / (x^2+1) , ## 2 terms for delta(x)
"delta3" = {tx2 <- (2/x)^2 ; -(1 - tx2*(1 - tx2)) / (x^2+1) }, ## 3 terms
stop("invalid 'method': ", method))
}
(S53meth <- eval(formals(eps_S53)$method))
## "direct" "log1pExct" "deltaFull" "delta1" "delta2" "delta3"
epsAllS53 <- function(x, meth = eval(formals(eps_S53)$method)) {
stopifnot(is.character(meth), is.numeric(x) || is(x,"mpfr"))
Rl <- lapply(meth, function(m) eps_S53(x, method=m))
new("mpfrMatrix", unlist(Rl),
Dim = c(length(x), length(meth)),
Dimnames = list(NULL, meth))
}
EA <- epsAllS53(x)
asNumeric(EA) # looks +- ok
## testing cbind() i.e., the Rmpfr-method a bit here !
chkDn <- function(r) stopifnot(identical(
dimnames(r),
list(NULL, c("x", "log2x", eval(formals(eps_S53)$method)))))
chkDn2 <- function(r) stopifnot(identical(
colnames(r)[-(1:2)], eval(formals(eps_S53)$method)))
chkDn ( r1 <- cbind(x, log2x = log2(x), EA))
chkDn2( r2 <- selectMethod("cbind","mpfr")(x, log2x = log2(x), EA) )
## are the same if names are *made* to match:
stopifnot(identical(r1, local({t <- r2; dimnames(t) <- dimnames(r1); t})))
chkDn2( r3 <- selectMethod("cbind","mpfr")(x, log2(x), epsAllS53(x)) )
chkDn2( r4 <- cbind(x, log2(x), epsAllS53(x)) )
stopifnot(identical(r1, local({t <- r4; dimnames(t) <- dimnames(r1); t})))
U_S53mat <- function(x) cbind(x, log2x = log2(x), epsAllS53(x))
chkDn(Rx <- U_S53mat(x))
dimnames(Rx) # .. all empty, i.e., list(NULL, NULL) ... even though epsAll..() gives nice colnames
##
if(FALSE) ## no longer needed !
## A version with correct colnames :
U_S53mat <- function(x) {
eA <- epsAllS53(x)
r <- cbind(x, log2(x), eA)
colnames(r) <- c("x", "log2.x", colnames(eA))
r
}
stopifnot(Rx[,-(1:2)] < 0)
op <- options(width = 100)
print(asNumeric(Rx), digits=5) ; options(op)
matplot(Rx[,"x"], - Rx[,-(1:2)], type = "b", log="xy")
## not useful really
opa <- par(ask=TRUE)
for(j in 3:ncol(Rx)) {
D <- Rx[,-c(1:2,j)] - Rx[,j]
neg <- mean(D < 0) > 1/2
matplot(Rx[,"x"], if(neg) -D else D,
main = sprintf("difference %s(Rx[,*] - Rx[, \"%s\"])",
if(neg) " -" else "", colnames(Rx)[j]),
type = "b", log="xy")
}
par(opa)
## rather:
str(epsMn <- rowMeans(Rx[,-(1:2)]))
matplot(Rx[,"x"], Rx[,-(1:2)] - epsMn, type = "b", log="x")
matplot(Rx[,"x"], abs(Rx[,-(1:2)] - epsMn), type = "b", log="xy") # don't see really
print(asNumeric(cbind(Rx[, 1:2], Rx[,-(1:2)] - epsMn)), digits=4)
if(FALSE) ## not used here anymore :
Mx <- cbind(x, log2x = log2(x), hSqrt = .5*(1-sqrt(1 + 8/x^2)), T1 = -2*x^-2, D = -2*x^-2 - .5*(1-sqrt(1 + 8/x^2)), D2= -2*x^-2*(1 - 2*x^-2) - .5*(1-sqrt(1 + 8/x^2)))
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.