R Under development (unstable) (2023-06-09 r84528 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 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. > library(robustbase) > > source(system.file("xtraR/styleData.R", package = "robustbase")) # -> smallD list of small datasets > str(smallD,, 20) List of 23 $ x0 : num(0) $ x1I : num Inf $ x1 : num 3 $ xII : num [1:2] -Inf Inf $ x2I : num [1:2] -Inf 9 $ x2 : int [1:2] 1 2 $ x3.2I: num [1:3] -Inf 9 Inf $ x3I : num [1:3] -Inf 9 11 $ x3 : num [1:3] 1 2 10 $ y : int [1:10] 1 2 3 4 5 6 7 8 9 10 $ y1 : num [1:11] 1 2 3 4 5 6 7 8 9 10 100 $ y. : int [1:11] 1 2 3 4 5 6 7 8 9 10 11 $ xC. : num [1:11] 1 1 1 1 1 1 1 1 1 1 1 $ yI : num [1:12] 1 2 3 4 5 6 7 8 9 10 100 Inf $ y2 : num [1:12] 1 2 3 4 5 6 7 8 9 10 100 1000 $ y1. : num [1:12] 1 2 3 4 5 6 7 8 9 10 11 100 $ xC1. : num [1:12] 1 1 1 1 1 1 1 1 1 1 1 10 $ xC : num [1:12] 1 1 1 1 1 1 1 1 1 1 1 1 $ yI. : num [1:13] 1 2 3 4 5 6 7 8 9 10 11 100 Inf $ y2. : num [1:13] 1 2 3 4 5 6 7 8 9 10 11 100 1000 $ xC2. : num [1:13] 1 1 1 1 1 1 1 1 1 1 1 10 100 $ xC1 : num [1:13] 1 1 1 1 1 1 1 1 1 1 1 1 10 $ xC2 : num [1:14] 1 1 1 1 1 1 1 1 1 1 1 1 10 100 > > lx <- lapply(smallD, + function(x) { + m <- mad(x) + hx <- + if(!is.na(m) && m > 0 && m != Inf) # in all these cases, MASS::huber() fails + MASS::huber(x) + else list(m=NA, s=NA) + hMx <- huberM(x) + list(loc = + c(median = median(x), + huber = hx$m, + huberM = hMx$m), + scale= + c(mad = m, + huber = hx$s, + huberM = hMx$s)) + }) > > > r <- list(mu = sapply(lx, function(x) x$loc), + s = sapply(lx, function(x) x$scale)) > r $mu x0 x1I x1 xII x2I x2 x3.2I x3I x3 y y1 y. xC. median NA Inf 3 NaN -Inf 1.5 9 9.000000 2.000000 5.5 6.000000 6 1 huber NA NA NA NA NA 1.5 NA 7.776102 2.611949 5.5 6.167169 6 NA huberM NA Inf 3 NaN -Inf 1.5 9 7.776102 2.611949 5.5 6.167169 6 1 yI y2 y1. xC1. xC yI. y2. xC2. xC1 xC2 median 6.500000 6.500000 6.500000 1 1 7.000000 7.000000 1 1 1 huber 6.834339 6.834339 6.606518 NA NA 7.213034 7.213034 NA NA NA huberM 6.834339 6.834339 6.606518 1 1 7.213034 7.213034 1 1 1 $s x0 x1I x1 xII x2I x2 x3.2I x3I x3 y y1 y. xC. mad NA NA 0 NA NA 0.7413 Inf 2.9652 1.4826 3.7065 4.4478 4.4478 0 huber NA NA NA NA NA 0.7413 NA 2.9652 1.4826 3.7065 4.4478 4.4478 NA huberM NA 0 0 0 0 0.7413 Inf 2.9652 1.4826 3.7065 4.4478 4.4478 0 yI y2 y1. xC1. xC yI. y2. xC2. xC1 xC2 mad 4.4478 4.4478 4.4478 0 0 4.4478 4.4478 0 0 0 huber 4.4478 4.4478 4.4478 NA NA 4.4478 4.4478 NA NA NA huberM 4.4478 4.4478 4.4478 0 0 4.4478 4.4478 0 0 0 > > cat('Time elapsed: ', proc.time(),'\n') # for ``statistical reasons'' Time elapsed: 0.25 0.06 0.25 NA NA > > proc.time() user system elapsed 0.25 0.06 0.25