R Under development (unstable) (2024-06-18 r86781 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 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.

> ### Have had cases where differences between large numbers lose precision, or even give Inf,
> ### which lead to NA
> 
> require(robustbase)
Loading required package: robustbase
> 
> source(system.file("xtraR/styleData.R", package = "robustbase"))# -> smallD, mkMx()
> 
> stopifnot(exprs = {
+     all.equal(scaleTau2(c(-4:4, 10000), consistency=FALSE),
+              (scaleTau2(c(-4:4, 1e300), consistency=FALSE) -> sT), # <- gave NaN, now fine !
+              tol = 1e-15) # even 0 (exact equality; Linux 64b)
+     all.equal(3.41103800034854, sT, tol = 1e-14) # seen 6.5e-16
+ })
> 
> exL <- c(
+     list( xNA  = c(NA, 1:6)
+        , xMe9 = c(-4:4, 1e9)
+        , xM   = c(-4:4, .Machine$double.xmax)
+        , xI   = c(-4:4, Inf)
+        , IxI  = c(-Inf, -4:4, Inf)
+        , IxI2 = c(-Inf, -4:4, Inf,Inf))
+     ,
+     mkMx(c(1e6, 1e9,
+            1e12, 1e14, 1e16,
+            1e20, 1e40, .Machine$double.xmax, Inf))
+ )
> 
> madL <- vapply(exL, mad, pi)
> ## Initially, scaleTau2() "works" but gives  NaN  everywhere -- now fine!
> sT2.L    <- vapply(exL, scaleTau2, FUN.VALUE=1, sigma0 = 1, consistency=FALSE)
> sT2.i2.L <- vapply(exL, scaleTau2, FUN.VALUE=1, sigma0 = 1, consistency=FALSE, iter = 2)
> sT2.i5.L <- vapply(exL, scaleTau2, FUN.VALUE=1, sigma0 = 1, consistency=FALSE, iter = 5)
> 
> cbind(madL, sT2.L)
                      madL    sT2.L
xNA                     NA       NA
xMe9          3.706500e+00 2.313007
xM            3.706500e+00 2.313007
xI            3.706500e+00 2.313007
IxI           4.447800e+00 2.412091
IxI2          4.447800e+00 2.440970
M1e06_nm1     1.482590e+06 2.664402
M1e06_neq     7.412993e+05 2.679394
M1e06_np1     1.037820e+01 2.706130
M1e09_nm1     1.482600e+09 2.664402
M1e09_neq     7.413000e+08 2.679394
M1e09_np1     1.037820e+01 2.706130
M1e12_nm1     1.482600e+12 2.664402
M1e12_neq     7.413000e+11 2.679394
M1e12_np1     1.037820e+01 2.706130
M1e14_nm1     1.482600e+14 2.664402
M1e14_neq     7.413000e+13 2.679394
M1e14_np1     1.037820e+01 2.706130
M1e16_nm1     1.482600e+16 2.664402
M1e16_neq     7.413000e+15 2.679394
M1e16_np1     1.037820e+01 2.706130
M1e20_nm1     1.482600e+20 2.664402
M1e20_neq     7.413000e+19 2.679394
M1e20_np1     1.037820e+01 2.706130
M1e40_nm1     1.482600e+40 2.664402
M1e40_neq     7.413000e+39 2.679394
M1e40_np1     1.037820e+01 2.706130
M1.8e308_nm1           Inf 2.664402
M1.8e308_neq 1.332630e+308 2.679394
M1.8e308_np1  1.037820e+01 2.706130
MInf_nm1               Inf 2.664402
MInf_neq               Inf 2.679394
MInf_np1      1.037820e+01 2.706130
> 
> stopifnot(exprs = {
+     is.na(madL [1])
+     is.na(sT2.L[1])
+     2.3 < sT2.L[-1]
+           sT2.L[-1] < 2.71
+ })
> 
> 
> xI <- exL$xI
> stopifnot(exprs = {
+     mad(exL$xI, constant = 1) == 2.5
+ })
> 
> stopifnot(exprs = {
+     all.equal(3.5471741782, scaleTau2(xI)) # gave NaN
+     ## These even gave  Error in ..... : NA/NaN/Inf in foreign function call (arg 1)
+     all.equal(3.5778,       Sn(xI))
+     all.equal(3.1961829592, Qn(xI))
+ })
> 
> ## From  example(mc)  {by MM} :
> 
> ## Susceptibility of the current algorithm to large outliers :
> dX10 <- function(X) c(1:5,7,10,15,25, X) # generate skewed size-10 with 'X'
> (Xs <- c(10,20,30, 60, 10^(2:10), 1000^(4:19), 1e6^c(10:20,10*(3:5)), Inf))
 [1]  1e+01  2e+01  3e+01  6e+01  1e+02  1e+03  1e+04  1e+05  1e+06  1e+07
[11]  1e+08  1e+09  1e+10  1e+12  1e+15  1e+18  1e+21  1e+24  1e+27  1e+30
[21]  1e+33  1e+36  1e+39  1e+42  1e+45  1e+48  1e+51  1e+54  1e+57  1e+60
[31]  1e+66  1e+72  1e+78  1e+84  1e+90  1e+96 1e+102 1e+108 1e+114 1e+120
[41] 1e+180 1e+240 1e+300    Inf
> 
> mc10x <- vapply(Xs, function(x) mc(dX10(x)), 1)
The default of 'doScale' is FALSE now for stability;
  set options(mc_doScale_quiet=TRUE) to suppress this (once per session) message
> ## now fixed:
> stopifnot(all.equal(c(4,6, rep(7,42))/12, mc10x))
> plot(Xs, mc10x, type="b", main = "mc( c(1:5,7,10,15,25, X) )", xlab="X", log="x")
> 
> ## so, Inf does work, indeed for mc()
> mcOld <- function(x, ..., doScale=TRUE) mc(x, doScale=doScale, c.huberize=Inf, ...)
> (x10I <- dX10(Inf))
 [1]   1   2   3   4   5   7  10  15  25 Inf
> set.seed(2020-12-04)# rlnorm(.)
> summary(xlN <- rlnorm(100))
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
 0.08955  0.55831  1.18474  1.91846  2.36054 14.88402 
> xII <- c(-Inf, xlN, Inf)
> stopifnot(exprs = {
+     all.equal(0.5,  print(mcOld(x10I)))
+     all.equal(7/12, print(mc   (x10I, doScale=TRUE ))) # was 0.5 before huberization
+     all.equal(7/12, print(mc   (x10I, doScale=FALSE)))
+     mcOld(xII) == 0
+     all.equal(0.3646680319, mc(xII))
+ })
[1] 0.5
[1] 0.5833333
[1] 0.5833333
> 
> proc.time()
   user  system elapsed 
   0.29    0.07    0.37