The Brobdingnag R package provides functions to enable the work with very large numbers. Most integration functions fail if large numbers temporarily occur. The R code here implements functions of the R package Brobdingnag to enable integration functions from R packages pracma and Bolstad2 to work with very large numbers. This is helpful if one uses functions like Gamma or others within integrals that easily reach values beyond normal tolerance ie. near zero or up to infinity. R's tolerance on a machine is
> 10^(307:309)
[1] 1e+307 1e+308 Inf
> 10^(-323:-324)
[1] 9.881313e-324 0.000000e+00
or more precisely
> .Machine$double.xmin
[1] 2.225074e-308
> .Machine$double.xmax
[1] 1.797693e+308
Another option is to transform all functions for which numerical integration is required to the log scale and then use a numerical integration method that is written for functions with log output. Normally this speeds up the process drastically compared to the usage of Brobdingnag while attaining the same accuracy, because the underlying R code is more or less exactly the same. However, it is slower than using non-log methods like integrate. Concretely, this means
fx <- function(x) sin(x)/cos(x)
becomes
fx.log <- function(x) log(sin(x)) - log(cos(x))
This conversion is no automatic process and for complicated integrals it requires a certain amount of effort and testing, esp. to make it efficient.
The third option is to use the package Rmpfr for very large numbers. It has a command called integrateR that works pretty fast.
| File | Description |
|---|---|
numinteg_accuracy.r |
check accuracy + worked example |
numinteg_brob.r |
numerical integration methods for brob objects (from Bolstad2: sintegral, from pracma: romberg, cotes, kronrod, quadgk, simpadpt, quadv, quadinf, quad, trapz, quadgr) |
numinteg_brob_calls.r |
example calls |
numinteg_log.r |
numerical integration methods for log functions (from cmna: trapez, romberg) |
numinteg_log_calls.r |
example calls |
simpsonrule_nlb.r |
general simpson rule for normal, log, and brob objects |
simpsonrule_nlb_calls.r |
example calls |
Original names just got a *.brob added at the end of the filename.
| Function | Description |
|---|---|
numinteg_accuracy.r |
Check accuracy + worked example |
ck.accuracy |
check accuracy |
numinteg_helper-funcs.r |
Helper functions for brob objects and log calculations |
list2vec.brob |
convert a brob list to a brob vector |
scalarprod.brob |
replace %*% scalarproduct that does not work for brob objects |
.llog2sub.short |
log(x-y) |
.llog.2add.shor |
log(x+y) |
numinteg_brob.r |
Numerical integration methods for brob objects |
sintegral.brob.parallel |
taken from Bolstad2 |
| taken from pracma: | |
romberg.brob |
Romberg |
cotes.brob |
Cotes |
integral.brob |
Kronrod, Simpson - but not Clenshaw |
quadgk.brob |
adaptive Gauss-Kronrod quadrature |
simpadpt.brob |
Simpson |
quadinf.brob |
infinite integrals |
quad.brob |
adaptive Simpson quadrature |
trapz.brob |
Trapez |
quadgr.brob |
Gaussian quadrature with Richardson extrapolation |
numinteg_log.r |
Numerical integration methods for log functions: |
| taken from cmna: | |
trap.nl |
Trapez |
romberg.nl |
Romberg |
simpsonrule_nlb.r |
Simpson rule |
simpsonrule.nlb |
Simpson for normal, log functions, and brob objects |
Note: Not all integration methods of pracma are covered. What is missing are:
- Clenshaw from
integral integral2integral3
The same is true for cmna.
Functions are called identical to the original version. Use as normal but add .brob to the name of the function call or as outlined above, and use a function that gives out a Brobdingnag object that acts as input for the integration function. For the log version you need a function that gives out a log value.
source("numinteg_helper-funcs.r")
source("numinteg_brob.r")
e.g.
fx <- function(x) sin(x)
f <- function(x) as.brob(sin(x))
lower <- 0
upper <- pi
maxit <- 25
tol <- 1e-15
romberg.brob(f, lower, upper, maxit=maxit, tol=tol)
as.brob(romberg(f=fx, a, b, maxit=maxit, tol=tol)$value)
The *_calls.r files contain examples for each function.
There is no special treatment to ensure accuracy embedded beyond what is already in the role models from pracma, cmna, and Bolstad2. If one is interested in that, you have to write your own tolerance error function based on known error functions of numerical integration methods. Another option and more pragmatic is to use functions for which the integration result can be obtained algebraically ie. via an analytical solution beyond any doubt. This can be compared to the output of the various nethods. Beyond the numerical precision of the machine one should not trust numbers anyway unless accuracy is a dedicated feature of the software and emulates accuracy. The file numinteg_accuracy.r contains a simple worked example.
# function to integrate algebraically and numerically
f <- function(x) sin(x)
# integral analytically solved
# sin(x) dx -> -cos(x) + constant
# integration boundaries
lower <- 0
upper <- pi
f.integ <- function(lower,upper) -cos(upper) -(-cos(lower))
# real value
int.real.v <- f.integ(lower=lower, upper=upper)
# numerical integration using method on log scale
f.log <- function(x) log(sin(x))
simpson.l.v <- exp(simpsonrule.nlb(f.log, lower, upper, method="log", Nsteps=Nsteps))
# check accuracy
ck.accuracy(comp=int.real.v, numinteg=simpson.l.v, methods=c("real value","simpson log"))
This results in:
Checking Numerical Integration Methods
--------------------------------------------
comparative method = real value
num. int. method = simpson log
comparative value = 2
num. integ. value = 2.0000000000000000
diff comp-num.integ. = 0.0000000000000000
comp. > num.integ. = FALSE
abs(diff) < tol = 1e-08 = TRUE
factor comp/num.integ. = 1.0000000000000000
---
Note: 'comparative' value can be a real value or from a different method
And now you can compare output results to a definite value to measure accuracy within the machine tolerance boundaries. Or you can compare two numerical integration methods. If you don't know whether your function can be solved algebraically, you can try with the online integration solver of wolframalpha.
For benchmarking, the R package microbenchmark is a good choice and easy to use. See the manpage of how to use it and some example calls here. As a short example we use the code from the accuracy check above:
> # benchmarking
> library(microbenchmark)
> Nsteps <- 1e3
> numinteg.benchm <- microbenchmark(
+ integrate(f, lower, upper)$v,
+ Rmpfr:::integrateR(f.Rmpfr, lower, upper)$value,
+ pracma:::trapz(x=sek, y=f(sek)),
+ simpsonrule.nlb(f, lower, upper, method="normal", Nsteps=Nsteps),
+ exp(simpsonrule.nlb(f.log, lower, upper, method="log", Nsteps=Nsteps)),
+ as.numeric( trapz.brob(x=sek, y=f.sek.brob.c) )
+ )
> numinteg.benchm
Unit: microseconds
expr min lq mean median uq max neval cld
integrate(f, lower, upper)$v 17.950 26.165 42.94091 46.615 55.640 71.14 100 a
Rmpfr:::integrateR(f.Rmpfr, lower, upper)$value 25575.446 26318.942 28044.54153 26948.197 28003.902 39713.58 100 b
pracma:::trapz(x = sek, y = f(sek)) 421.650 436.935 465.50396 445.700 454.520 829.30 100 a
simpsonrule.nlb(f, lower, upper, method = "normal", Nsteps = Nsteps) 62.140 71.510 90.17502 91.560 98.245 269.67 100 a
exp(simpsonrule.nlb(f.log, lower, upper, method = "log", Nsteps = Nsteps)) 1403.640 1477.560 1687.81533 1521.095 1572.291 15313.48 100 c
as.numeric(trapz.brob(x = sek, y = f.sek.brob.c)) 7607.822 7882.082 8339.16269 8047.017 8205.902 21620.37 100 d
Note:
Different methods are mixed here (e.g. Simpson rule, Trapez, Romberg (integrateR), adaptive quadrature (integrate)). Additionally, the number of steps Nsteps differs between simpsonrule.nlb, trapz, integrate, and integrateR esp. if a cut-off criterium is chosen based on some tolerance value. Thus, causes of any differences in computer time are difficult to identify precisely. For serious comparison one should standardize all influences so that a difference is indeed caused by either the input method (like log)() or as.brob()/ brob()) and not by the algorithm itself (like Simpson, Trapez, Romberg, etc.) or just by the number of steps Nsteps used to calculate the integral. This is valid for all other comparisons. As an ideal all influences should (must) be held constant except for what you want to compare. Beyond that, benchmarking is not everything -- accuracy must be considered too if you choose a method. Both together allow for a compromise if you work with very large numbers.
Thus, we can repeat the benchmarking above but we will use only one method simpsonrule.nlb but with different input: normal, log, brob. First, we check for accuracy:
> Nsteps <- 1e3
> # repeat for one method
> # add an alternative function
> f.brob2 <- function(x) brob(log(sin(x)))
>
> simpsonrule.nlb.n <- simpsonrule.nlb(f, lower, upper, method="normal", Nsteps=Nsteps)
> simpsonrule.nlb.l <- exp(simpsonrule.nlb(f.log, lower, upper, method="log", Nsteps=Nsteps))
> simpsonrule.nlb.b <- as.numeric(simpsonrule.nlb(f.brob, lower, upper, method="brob", Nsteps=Nsteps))
> simpsonrule.nlb.b2 <- as.numeric(simpsonrule.nlb(f.brob2, lower, upper, method="brob", Nsteps=Nsteps))
>
> ck.accuracy(comp=int.real.v, numinteg=simpsonrule.nlb.n, methods=c("real","simpsonrule.nlb - normal"))
Checking Numerical Integration Methods
--------------------------------------------
comparative method = real
num. int. method = simpsonrule.nlb - normal
comparative value = 2
num. integ. value = 1.9999991775331358
diff comp-num.integ. = 0.0000008224668642
comp. > num.integ. = TRUE
abs(diff) < tol = 1e-08 = FALSE
factor comp/num.integ. = 1.0000004112336012
---
Note: 'comparative' value can be a real value or from a different method
> ck.accuracy(comp=int.real.v, numinteg=simpsonrule.nlb.l, methods=c("real","simpsonrule.nlb - log"))
Checking Numerical Integration Methods
--------------------------------------------
comparative method = real
num. int. method = simpsonrule.nlb - log
comparative value = 2
num. integ. value = 2.0000000000000675
diff comp-num.integ. = -0.0000000000000675
comp. > num.integ. = FALSE
abs(diff) < tol = 1e-08 = TRUE
factor comp/num.integ. = 0.9999999999999662
---
Note: 'comparative' value can be a real value or from a different method
> ck.accuracy(comp=int.real.v, numinteg=simpsonrule.nlb.b, methods=c("real","simpsonrule.nlb - brob"))
Checking Numerical Integration Methods
--------------------------------------------
comparative method = real
num. int. method = simpsonrule.nlb - brob
comparative value = 2
num. integ. value = 2.0000000000000657
diff comp-num.integ. = -0.0000000000000657
comp. > num.integ. = FALSE
abs(diff) < tol = 1e-08 = TRUE
factor comp/num.integ. = 0.9999999999999671
---
Note: 'comparative' value can be a real value or from a different method
> ck.accuracy(comp=int.real.v, numinteg=simpsonrule.nlb.b2, methods=c("real","simpsonrule.nlb - brob2"))
Checking Numerical Integration Methods
--------------------------------------------
comparative method = real
num. int. method = simpsonrule.nlb - brob2
comparative value = 2
num. integ. value = 2.0000000000000657
diff comp-num.integ. = -0.0000000000000657
comp. > num.integ. = FALSE
abs(diff) < tol = 1e-08 = TRUE
factor comp/num.integ. = 0.9999999999999671
---
Note: 'comparative' value can be a real value or from a different method
and then we compare the speed:
> # just compare within one method but with different input
> numinteg.benchm2 <- microbenchmark(
+ simpsonrule.nlb(f, lower, upper, method="normal", Nsteps=Nsteps),
+ simpsonrule.nlb(f.log, lower, upper, method="log", Nsteps=Nsteps),
+ simpsonrule.nlb(f.brob, lower, upper, method="brob", Nsteps=Nsteps),
+ simpsonrule.nlb(f.brob2, lower, upper, method="brob", Nsteps=Nsteps)
+ )
> numinteg.benchm2
Unit: microseconds
expr min lq mean median uq max neval cld
simpsonrule.nlb(f, lower, upper, method = "normal", Nsteps = Nsteps) 61.64 66.79 81.06962 78.810 92.600 112.64 100 a
simpsonrule.nlb(f.log, lower, upper, method = "log", Nsteps = Nsteps) 1398.98 1485.38 1693.13454 1525.745 1609.541 12096.86 100 a
simpsonrule.nlb(f.brob, lower, upper, method = "brob", Nsteps = Nsteps) 231282.95 240144.51 252877.14441 251100.391 258555.633 515170.27 100 b
simpsonrule.nlb(f.brob2, lower, upper, method = "brob", Nsteps = Nsteps) 157890.80 167820.18 176322.08066 175793.131 183799.018 211417.95 100 c
Note:
As we can see the latter two ie. f.brob and f.brob2 lead to exactly the same results, use the same calculation procedure but differ greatly in time. The only difference is as.brob(x) vs. brob(log(x)) whereas the first one is much slower than the second one.
Using parallel for multi-threading does not lead to expected results -- see for yourself. It may be helpful if a function to be integrated numerically is highly complicated and it requires a lot of repetiions. Otherwise, trials showed no advantage of using parallel over normal computing. Rather, the opposite was the case: The initial overhead of parallel computing required more time than the normal time required to compute everything. If this is repeated, the same story starts again. So, there is no gain at all at this stage. Future efforts will show whether this can be changed.
Cover all integration methods and add them as multi-threaded nevertheless to find out how to integrate parallel computing better and more efficiently to speed up the process.
See license file included - in short it is GPL >=2 - in accordance to the packages used and the license associated with them.
All R scripts should work under R >=v3.
The R code was tested carefully and cross-checked against various public available results and manpages to ensure proper results. The example calls can be used to test for accuracy. However, it is provided "as is". Use common sense to compare results with expectations. NO WARRANTY of any kind is involved here. There is no guarantee that the software is free of error or consistent with any standards or even meets your requirements. Do not use the software or rely on it to solve problems if incorrect results may lead to hurting or injurying living beings of any kind or if it can lead to loss of property or any other possible damage to the world. If you use the software in such a manner, you are on your own and it is your own risk.