Nothing
#### Tests for "Functionals": Root finding, Optimization, Integration, etc
stopifnot(require("Rmpfr"))
(f.chk <- system.file("check-tools.R", package="Rmpfr", mustWork=TRUE))
source(f.chk, keep.source=FALSE)
## -> assert.EQ.() showSys.time()
options(warn = 1)# warnings *immediately*
(doExtras <- Rmpfr:::doExtras())
### 1. Integration -----------------------------------------------
## Example from Lauren K, June 2014 (~/R/MM/Pkg-ex/Rmpfr/integrateR-LaurenK.R):
beta0 <- 0.05
beta1 <- 0.05
(tau <- sqrt(0.01*0.05*0.95/0.99))# = 0.0219..
##
Z00 <- 9
Z01 <- 1
Z10 <- 18
Z11 <- 2
N <- Z00+Z01+Z10+Z11
integrand <- function(u) {
ee.u <- exp(-u^2/2)/(sqrt(2*pi)*tau)
b0u <- beta0 + tau*u
b1u <- beta1 + b0u # == beta0+beta1+ tau*u
ee.u ^ (Z00+Z01 + Z10+Z11) *
(1-b0u)^Z00 * b0u ^ Z01 *
(1-b1u)^Z10 * b1u ^ Z11
}
## MM: note how the integrand() function looks:
op <- par(mfcol=c(3,1), mgp=c(1.5, .6, 0), mar=.1+c(3,3,0,0))
ax1 <- function(a,b) axis(1, at=c(a,b), col=2, col.axis=2, tcl= +3/4, mgp=c(3,0,0))
curve(integrand, -5,5, n=1000)
cc <- adjustcolor(2, 1/4) ## look closer:
ep <- .01; rect(-3, -ep, 0, +ep, col=cc, border=cc); ax1(-3,0)
curve(integrand, -3,0, n=1000, ylim = c(-ep,ep))
## but now look really closely :
ep <- .001; rect(-3, -ep, -2, +ep, col=cc); ax1(-3,-2)
curve(integrand, -3,-2, n=1000, ylim = c(-ep, ep))
par(op)
(I1 <- integrate(integrand,lower = -100, upper = 100))
(I1. <- integrate(integrand,lower = -100, upper = 100, rel.tol = 1e-14))
showSys.time(I2 <- integrateR(integrand, lower = -100, upper = 100))
I2 ## ... Warning ‘no convergence up to order 13’
## Solaris Sparc (2014-06, CRAN checks); thanks Brian: print(I2[1:2], digits=15)
I2.Solaris <- list(value = 1.3963550396006e+33, abs.error = 1.79487857486724e+28)
I.db <- list(value = 1.39635503960059e+33, abs.error = 1.79487857478077e+28)
stopifnot(
all.equal(I2[1:2], I.db, tol = 1e-10)# Solaris SPARC needs at least 4.8e-11
)
## Now using high accuracy
showSys.time(I3 <- integrateR(integrand, lower = mpfr(-100, precBits=256), upper = 100))
## much slower but not better (and not worse)
I3
assert.EQ.(sapply(I3[1:2], asNumeric), unlist(I.db))
## Really get better when decreasing the integration interval
## from [-100, 100] to [-10, 10] ... which should give "the same"
showSys.time(I4 <- integrateR(integrand, lower = mpfr(-10, precBits=256), upper = 10,
ord = 15, verbose=TRUE))
## ~ 6.6 sec [lynne 2013]
I4
## on the left side, there is "nothing" (and negative, as we know!):
showSys.time(I0.1 <- integrateR(integrand, lower = mpfr(-1000, precBits=256),
upper = -10, ord= 11, verbose=TRUE))
showSys.time(I0.2 <- integrateR(integrand, lower = mpfr(10, precBits=256),
upper = 1000, ord= 11, verbose=TRUE))
I0.1
I0.2
I4
## Integral [-1000, +1000 ] = Int[-1000, -10] + Int[-10, +10] + Int[+10, +1000]:
I4 $value + I0.1 $value + I0.2 $value
## but this is really the same as just the middle:
stopifnot(I4 $value + I0.1 $value + I0.2 $value
== I4 $value)
value <- I4$value; delta <- I4$abs.err
nDig <- -asNumeric(log10(delta/value))
cat("Correct number of digits: ", round(nDig, 2),"\n",
"Integral I = ", format(I4$value, digits = ceiling(nDig)),
" (last change change= ", format(delta, digits = 7),")\n",
"integrate(.) = ", format(I1 $value, digits = 22),"\n",
"integrate(., rtol=1e-15)= ", format(I1.$value, digits = 22),"\n", sep="")
### 2. Root Finding ----------------------------------------------
### 3. Optimization / Minimization, .. ---------------------------
cat('Time elapsed: ', proc.time(),'\n') # "stats"
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.