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