R Under development (unstable) (2026-09-25 r90590 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2026 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> #### Original in ~/R/Pkgs/DistributionUtils/inst/unitTests/runit.incompleteBesselK.R
> 
> ### "Translated" from RUnit unit tests to simple base R tests by Martin Maechler
> ##  checkEquals(....)              ==> stopifnot(all.equal(<....>))
> ##  checkTrue  (...., msg = <msg>) ==> stopifnot(<msg> = <....>)
> 
> ### Unit tests of function incompleteBesselK and incompleteBesselKR
> 
> ## Author: David Scott, Date: 5 Jan 2012
> 
> ## Original version of the test file which was misplaced
> 
> library(Bessel)
> ## Data from Harris (2008)
> HarrisCase1 <- c(2.225310761266469,
+                  0.213894166822940,
+                  0.054503469799701,
+                  0.023253121507708,
+                  0.013042750996080,
+                  0.008567534990649,
+                  0.006208676806601,
+                  0.004801085238177,
+                  0.003884072049627,
+                  0.003246798003149)
> 
> HarrisCase2 <- 0.000012249987981
> 
> HarrisCase3 <- 0.00000041500106421228
> 
> HarrisCase4 <- 0.000528504325244
> 
> ## Calculations for Harris Case 1
> numIBF <- numeric(10)
> for (i in 0:9){
+     ibf <- incompleteBesselK(0.01, 4, i, nmax = 100)
+     numIBF[i + 1] <- ibf
+ }
> maxDiff <- max(abs(numIBF - HarrisCase1))
> ## stopifnot(paste("Harris Case 1: maxDiff =", maxDiff) = maxDiff <= 10^(-13))
> if(!(maxDiff <= 10^(-13))) stop(paste("Harris Case 1: maxDiff =", maxDiff))
> 
> ## Calculations for Harris Case 2
> ibf <- incompleteBesselK(4.95, 5, 2)
> maxDiff <- max(abs(ibf - HarrisCase2))
> stopifnot("Harris Case 2" = maxDiff <= 10^(-13))
> 
> ## Calculations for Harris Case 3
> ibf <- incompleteBesselK(10, 2, 6)
> maxDiff <- max(abs(ibf - HarrisCase3))
> stopifnot("Harris Case 3" = maxDiff <= 10^(-11))
> 
> ## Calculations for Harris Case 4
> ibf <- incompleteBesselK(3.1, 2.6, 5)
> maxDiff <- max(abs(ibf - HarrisCase4))
> ## stopifnot(paste("Harris Case 4: maxDiff =", maxDiff) = maxDiff <= 10^(-13))
> if(!(maxDiff <= 10^(-13))) stop(paste("Harris Case 4: maxDiff =", maxDiff))
> 
> ## Newer version of the test file: not necessarily better
> ## But includes testing of the pure R version
> 
> ## Values given by Harris (2008)
> ## MM: should ask for more precision than default sqrt(.Machine$double.eps) = 1.4901e-8
> stopifnot(all.equal(incompleteBesselK (0.01, 4, 0), 2.225310761266469))
> stopifnot(all.equal(incompleteBesselKR(0.01, 4, 0), 2.225310761266469))
> 
> stopifnot(all.equal(incompleteBesselK (0.01, 4, 9), 0.003246798003149))
> stopifnot(all.equal(incompleteBesselKR(0.01, 4, 9), 0.003246798003149))
> 
> stopifnot(all.equal(incompleteBesselK (4.95, 5, 2), 0.000012249987981))
> stopifnot(all.equal(incompleteBesselKR(4.95, 5, 2), 0.000012249987981))
> 
> stopifnot(all.equal(incompleteBesselK (10, 2, 6), 0.0000004150045864731308))
> stopifnot(all.equal(incompleteBesselKR(10, 2, 6), 0.0000004150045864731308))
> 
> stopifnot(all.equal(incompleteBesselK (3.1, 2.6, 5), 0.000528504325244))
> stopifnot(all.equal(incompleteBesselKR(3.1, 2.6, 5), 0.000528504325244))
> 
> ### Check values when x > y using numeric integration:
> integrBess <- function(t, x, y, nu)  t^-(nu+1) * exp(-x*t - y/t)
> 
> nus <- seq(0, 12, by = 1/4)
> x. <- 4
> y. <- 0.01
> 
> numIBF <- sapply(nus, incompleteBesselK, x = x., y = y.)
> integIBF <- sapply(nus, integrate, f = integrBess, lower = 1, upper = Inf, x = x., y = y.)
> ## --> a  5 x #{nus}  matrix of lists
> intIBF <- as.numeric(integIBF["value", ])
> errIBF <- as.numeric(integIBF["abs.error", ])
> numIBF - intIBF # -1.257989e-11 -1.154437e-11 -9.600298e-12 -7.442554e-12 -5.422429e-12 .....
 [1] -1.257989e-11 -1.154437e-11 -9.600298e-12 -7.442554e-12 -5.422429e-12
 [6] -3.746299e-12 -2.440144e-12 -1.480996e-12 -8.149757e-13 -3.791242e-13
[11] -1.135398e-13  3.294543e-14  1.008698e-13  1.204488e-13  1.130302e-13
[16]  9.278125e-14  6.856234e-14  4.538253e-14  2.567955e-14  1.029515e-14
[21] -8.554355e-16 -8.347706e-15 -1.290613e-14 -1.528725e-14 -1.614052e-14
[26] -1.599220e-14 -1.520659e-14 -3.181223e-14 -2.920082e-14 -2.665619e-14
[31] -2.422606e-14 -2.197244e-14 -1.986627e-14 -1.797520e-14 -1.623333e-14
[36] -1.468205e-14 -1.327172e-14 -1.202445e-14 -2.860104e-14 -2.593000e-14
[41] -2.355126e-14 -2.135488e-14 -1.933978e-14 -1.752309e-14 -1.587380e-14
[46] -1.434096e-14 -1.298853e-14 -1.177400e-14 -1.061369e-14
> stopifnot(identical(length(nus), as.integer(table(unlist(integIBF["message",])))),
+           max(abs(numIBF - intIBF)) < 1e-10,
+           abs(numIBF - intIBF) < errIBF)
> 
> plot(nus, numIBF, type = "o", col = 2, cex = 1/2)
> 
> 
> proc.time()
   user  system elapsed 
   0.57    0.09    0.65