1### Have had cases where differences between large numbers lose precision, or even give Inf,
2### which lead to NA
3
4require(robustbase)
5
6stopifnot(exprs = {
7    all.equal(scaleTau2(c(-4:4, 10000), consistency=FALSE),
8             (scaleTau2(c(-4:4, 1e300), consistency=FALSE) -> sT), # <- gave NaN, now fine !
9             tol = 1e-15) # even 0 (exact equality; Linux 64b)
10    all.equal(3.41103800034854, sT, tol = 1e-14) # seen 6.5e-16
11})
12
13mkMx <- function(M, ngood = 10, left = floor(ngood/3)) {
14    stopifnot(is.numeric(ngood), ngood >= 3,
15              is.numeric(M), length(M) == 1L, M >= 1000,
16              is.numeric(left), 0 <= left, left <= ngood)
17    right <- ngood-left
18    res <- list(
19        c(rep(-M, left), seq_len(ngood - 1L), rep(M, right)) # < 50% "good"
20      , c(rep(-M, left), seq_len(ngood     ), rep(M, right)) # half  "good"
21      , c(rep(-M, left), seq_len(ngood + 1L), rep(M, right)) # > 50% "good"
22    )
23    nM <- gsub("[-+]", "", formatC(M, digits=2, width=1))
24    names(res) <- paste0("M", nM,"_n", c("m1", "eq", "p1"))
25    res
26}
27
28exL <- c(
29    list( xNA  = c(NA, 1:6)
30       , xMe9 = c(-4:4, 1e9)
31       , xM   = c(-4:4, .Machine$double.xmax)
32       , xI   = c(-4:4, Inf)
33       , IxI  = c(-Inf, -4:4, Inf)
34       , IxI2 = c(-Inf, -4:4, Inf,Inf))
35  ##
36, mkMx(M = .Machine$double.xmax)
37, mkMx(M = 1e6)
38, mkMx(M = 1e9)
39, mkMx(M = 1e12)
40, mkMx(M = 1e14)
41, mkMx(M = 1e16)
42, mkMx(M = 1e20)
43, mkMx(M = 1e40)
44)
45
46madL <- vapply(exL, mad, pi)
47## Initially, scaleTau2() "works" but gives  NaN  everywhere -- now fine!
48sT2.L    <- vapply(exL, scaleTau2, FUN.VALUE=1, sigma0 = 1, consistency=FALSE)
49sT2.i2.L <- vapply(exL, scaleTau2, FUN.VALUE=1, sigma0 = 1, consistency=FALSE, iter = 2)
50sT2.i5.L <- vapply(exL, scaleTau2, FUN.VALUE=1, sigma0 = 1, consistency=FALSE, iter = 5)
51
52cbind(madL, sT2.L)
53
54stopifnot(exprs = {
55    is.na(madL [1])
56    is.na(sT2.L[1])
57    2.3 < sT2.L[-1]
58          sT2.L[-1] < 2.71
59})
60
61
62xI <- exL$xI
63stopifnot(exprs = {
64    mad(exL$xI, constant = 1) == 2.5
65})
66
67## FIXME: should not give NaN :
68scaleTau2(xI)
69
70## FIXME: even give  Error in ..... : NA/NaN/Inf in foreign function call (arg 1)
71try( Sn(xI) )
72try( Qn(xI) )
73
74## From  example(mc)  {by MM} :
75
76## Susceptibility of the current algorithm to large outliers :
77dX10 <- function(X) c(1:5,7,10,15,25, X) # generate skewed size-10 with 'X'
78(Xs <- c(10,20,30, 60, 10^(2:10), 1000^(4:19), 1e6^c(10:20,10*(3:5)), Inf))
79
80(mc10x <- vapply(Xs, function(x) mc(dX10(x)), 1))
81plot(Xs, mc10x, type="b", main = "mc( c(1:5,7,10,15,25, X) )", xlab="X", log="x")
82
83##--- FIXME: the above must change!
84
85## so, Inf does work, indeed for mc()
86dX10(Inf)
87set.seed(2020-12-04)
88stopifnot(exprs = {
89    is.finite(mc(dX10(Inf))) # 0.5 currently
90    mc(c(-Inf, rlnorm(100), Inf)) == 0
91})
92
93