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. > #### Test special cases for ltsReg() > > library(robustbase) > > ## Platform - and other such info -- so we find it in old saved outputs > .libPaths() [1] "D:/temp/RtmpeGLmbV/RLIBS_1b8a8f0e2c25" [2] "D:/RCompile/recent/R/library" > SysI <- Sys.info() > structure(Sys.info()[c(4,5,1:3)], class="simple.list") _ nodename CRANWIN3 machine x86-64 sysname Windows release Server x64 version build 20348 > library(lib.loc = .libPaths()[1]) # the "R CMD check specific tmp-library" > sessionInfo() R Under development (unstable) (2023-06-09 r84528 ucrt) Platform: x86_64-w64-mingw32/x64 Running under: Windows Server 2022 x64 (build 20348) Matrix products: default locale: [1] LC_COLLATE=C LC_CTYPE=German_Germany.utf8 [3] LC_MONETARY=C LC_NUMERIC=C [5] LC_TIME=C time zone: Europe/Berlin tzcode source: internal attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] robustbase_0.99-0 loaded via a namespace (and not attached): [1] DEoptimR_1.0-14 compiler_4.4.0 > c(robustbase = packageDescription("robustbase")$Built, + DEoptimR = packageDescription("DEoptimR")$Built) robustbase "R 4.4.0; x86_64-w64-mingw32; 2023-06-12 21:50:06 UTC; windows" DEoptimR "R 4.4.0; ; 2023-06-11 23:51:19 UTC; windows" > if(SysI[["sysname"]] == "Linux" && require("sfsmisc")) local({ + nn <- names(.Sc <- sfsmisc::Sys.cpuinfo()) + nn <- names(.Sc <- .Sc[nn != "flags"]) + print(.Sc[grep("\\.[0-9]$", nn, invert=TRUE)]) + }) > > ### 1) p = 1 ---------------------------------------------------- > set.seed(1) > x <- c(rnorm(50),100, 1e10) > (r1 <- ltsReg(x ~ 1)) # failed in Valentin's 1.0-3 (pre-version) Call: ltsReg.formula(formula = x ~ 1) Coefficients: [1] 0.1004 Scale estimate 0.9435 > summary(r1) Call: ltsReg.formula(formula = x ~ 1) Residuals (from reweighted LS): Min 1Q Median 3Q Max -2.0898 -0.3668 0.0000 0.6069 1.4948 Coefficients: Estimate Std. Error t value Pr(>|t|) [1,] 0.1004 0.1101 0.912 0.366 Residual standard error: 0.7707 on 48 degrees of freedom > (r1. <- ltsReg(y = x)) Call: ltsReg.default(y = x) Coefficients: [1] 0.1004 Scale estimate 0.9435 > i1 <- 15:17; ii <- (1:20)[-i1] > UN <- function(lis) lapply(lis, unname) > dimnames(r1.$X)[1] <- dimnames(r1$X)[1] > stopifnot(all.equal( r1[ii], r1.[ii], tolerance= 1e-15), + all.equal(UN(r1[i1]), UN(r1.[i1]), tolerance= 1e-15)) > > ## intercept=FALSE, p > 1 -- coefficients were switched once > n <- 100; theta <- c(x=10, x2=40) > set.seed(7) > X <- cbind(x = rt(n, 4), x2 = rnorm(n)) > dat <- data.frame(X, y = X %*% theta + rt(n, df=3)/10) > summary(M <- ltsReg(y ~ . -1, data = dat)) Call: ltsReg.formula(formula = y ~ . - 1, data = dat) Residuals (from reweighted LS): Min 1Q Median 3Q Max -0.287883 -0.047986 0.003272 0.083622 0.292822 Coefficients: Estimate Std. Error t value Pr(>|t|) x 9.991585 0.008511 1174 <2e-16 *** x2 40.011634 0.012009 3332 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.1167 on 90 degrees of freedom Multiple R-Squared: 1, Adjusted R-squared: 1 F-statistic: 6.911e+06 on 2 and 90 DF, p-value: < 2.2e-16 > stopifnot(all.equal(coef(M), theta, tolerance = 1e-3)) > > ## with alpha = 1 > (r1.1 <- ltsReg(x ~ 1, alpha = 1)) Call: ltsReg.formula(formula = x ~ 1, alpha = 1) Coefficients: [1] 2.059 Scale estimate 15.11 > summary(r1.1) Call: ltsReg.formula(formula = x ~ 1, alpha = 1) Residuals (from reweighted LS): Min 1Q Median 3Q Max -4.274 -2.387 -1.803 -1.295 0.000 Coefficients: Estimate Std. Error t value Pr(>|t|) [1,] 2.0593 0.3035 6.784 1.43e-08 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 2.146 on 49 degrees of freedom > > ### 1b) p = 1, constant scale > (rc <- ltsReg(y = rep(1,12))) Call: ltsReg.default(y = rep(1, 12)) Coefficients: [1] 1 Scale estimate 0 > str(rc) List of 22 $ method : chr "Univariate location and scale estimation.\nMore than half of the data are equal!" $ best : NULL $ coefficients : num 1 $ alpha : num 0.5 $ quan : num 7 $ raw.coefficients: num 1 $ raw.resid : num [1:12] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... $ raw.weights : num [1:12] 1 1 1 1 1 1 1 1 1 1 ... $ raw.scale : num 0 $ scale : num 0 $ crit : num 0 $ resid : num [1:12] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... $ rsquared : num 0 $ intercept : logi TRUE $ lts.wt : num [1:12] 1 1 1 1 1 1 1 1 1 1 ... $ residuals : num [1:12] 0 0 0 0 0 0 0 0 0 0 ... $ fitted.values : num [1:12] 1 1 1 1 1 1 1 1 1 1 ... $ Y : num [1:12] 1 1 1 1 1 1 1 1 1 1 ... $ X : num [1:12, 1] 1 1 1 1 1 1 1 1 1 1 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : NULL .. ..$ : NULL $ raw.cnp2 : num [1:2] 4.97 1.19 $ cnp2 : num [1:2] 1 1 $ call : language ltsReg.default(y = rep(1, 12)) - attr(*, "class")= chr "lts" > summary(rc) Call: ltsReg.default(y = rep(1, 12)) Residuals (from reweighted LS): Min 1Q Median 3Q Max 0 0 0 0 0 Coefficients: Estimate Std. Error t value Pr(>|t|) [1,] 1 0 Inf <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0 on 11 degrees of freedom > ## with alpha = 1 > (rc1 <- ltsReg(y = rep(1,12), alpha = 1)) Call: ltsReg.default(y = rep(1, 12), alpha = 1) Coefficients: [1] 1 Scale estimate 0 > summary(rc1) Call: ltsReg.default(y = rep(1, 12), alpha = 1) Residuals (from reweighted LS): Min 1Q Median 3Q Max 0 0 0 0 0 Coefficients: Estimate Std. Error t value Pr(>|t|) [1,] 1 0 Inf <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0 on 11 degrees of freedom > stopifnot(residuals(rc) == 0, all.equal(unname(coef(rc )), 1), + residuals(rc1) == 0, all.equal(unname(coef(rc1)), 1)) > > ### 2) alpha = 1 : classical estimates --- for general cases -------- > > > cat('Time elapsed: ', proc.time(),'\n') # for ``statistical reasons'' Time elapsed: 0.25 0.14 0.39 NA NA > > proc.time() user system elapsed 0.25 0.14 0.39