R Under development (unstable) (2026-02-18 r89435 ucrt) -- "Unsuffered Consequences" Copyright (C) 2026 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(robustX) > library(robustbase) > (newRB <- (packageVersion("robustbase") >= "0.99")) [1] TRUE > > sessionInfo() R Under development (unstable) (2026-02-18 r89435 ucrt) Platform: x86_64-w64-mingw32/x64 Running under: Windows Server 2022 x64 (build 20348) Matrix products: default LAPACK version 3.12.1 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-7 robustX_1.2-8 loaded via a namespace (and not attached): [1] DEoptimR_1.1-4 compiler_4.6.0 > packageDescription("robustX") Package: robustX Type: Package Title: 'eXtra' / 'eXperimental' Functionality for Robust Statistics Version: 1.2-8 VersionNote: Last CRAN: 1.2-7 on 2023-06-16 Date: 2026-02-20 Authors@R: c(person("Martin","Maechler", role=c("aut","cre"), email="maechler@stat.math.ethz.ch", comment = c(ORCID = "0000-0002-8685-9910")) , person("Werner A.", "Stahel", role="aut", email="stahel@stat.math.ethz.ch") , person("Rolf", "Turner", role="ctb", email="r.turner@auckland.ac.nz", comment = "reclas()") , person("Ueli", "Oetliker", role="ctb", comment = "original version of BACON() and mvBACON for S+") , person("Tobias", "Schoch", role="ctb", comment = "init.sel=\"V2\" for BACON; fix alpha") ) Maintainer: Martin Maechler URL: https://robustbase.R-forge.R-project.org/, https://R-forge.R-project.org/R/?group_id=59, https://R-forge.R-project.org/scm/viewvc.php/pkg/robustX/?root=robustbase, svn://svn.r-forge.r-project.org/svnroot/robustbase/pkg/robustX BugReports: https://R-forge.R-project.org/tracker/?atid=302&group_id=59 Description: Robustness -- 'eXperimental', 'eXtraneous', or 'eXtraordinary' Functionality for Robust Statistics. Hence methods which are not well established, often related to methods in package 'robustbase'. Amazingly, 'BACON()', originally by Billor, Hadi, and Velleman (2000) has become established in places. The "barrow wheel" `rbwheel()` is from Stahel and Mächler (2009) . Imports: grDevices, graphics, stats, utils, robustbase (>= 0.92-3) Suggests: MASS, lattice, pcaPP Enhances: ICS License: GPL (>= 2) Encoding: UTF-8 NeedsCompilation: no Packaged: 2026-02-21 22:05:11 UTC; maechler Author: Martin Maechler [aut, cre] (ORCID: ), Werner A. Stahel [aut], Rolf Turner [ctb] (reclas()), Ueli Oetliker [ctb] (original version of BACON() and mvBACON for S+), Tobias Schoch [ctb] (init.sel="V2" for BACON; fix alpha) Built: R 4.6.0; ; 2026-02-23 10:45:19 UTC; windows -- File: D:/RCompile/CRANincoming/R-devel/lib/robustX/Meta/package.rds > (ourBLAS <- grepl(print(normalizePath(R.home())), + normalizePath(extSoftVersion()[["BLAS"]]), fixed = TRUE)) [1] "D:\\RCompile\\recent\\R" [1] FALSE Warning message: In normalizePath(path.expand(path), winslash, mustWork) : path[1]="": Das System kann den angegebenen Pfad nicht finden > ## need extended precision (typically *includes* 64-bit): > doCheck <- (.Machine$sizeof.longdouble >= 16) > cat("doCheck (= have long double):", doCheck,"\n") doCheck (= have long double): TRUE > > if(!dev.interactive(orNone=TRUE)) pdf("cov-ex.pdf") > > ## report.and.stop.if.not.all.equal > report.stopifnot.all.eq <- function(a,b, tol, ...) { + call <- sys.call() + ae <- all.equal(a,b, tolerance=tol, ...) + call[[1]] <- quote(all.equal) + if(!isTRUE(ae)) + stop(sprintf("Not %s:\n%s\n\n", deparse(call), + paste(ae, collapse="\n")), + call.=FALSE) + ## else + TRUE + } > > UN <- function(L) lapply(L, unname) > > chk.NN.new.old <- function(cNew, cNold, tol = 2e-15, tol.1 = 20*tol) { + stopifnot(is.list(cNold$innc), length(n.i <- names(cNold$innc)) == 4) + cat("classification accordance matrix:\n") + print(table(new = cNew $classification, + old = cNold$classification)) + report.stopifnot.all.eq(UN(cNew [1:4]), + UN(cNold[1:4]), tol=tol.1) & + report.stopifnot.all.eq(cNew $innc[n.i], + cNold$innc[n.i], tol=tol) + } > > summ.NN <- function(cNN, digits = 3) { + cbind(class = cNN$classification, + pprob = round(cNN$postprob, digits), + incc.p= round(cNN$innc$postprob, digits)) + } > > > ## the original definition (2003) "plus" uniqEV > covNN.1 <- function(X) robustX:::covNNC1(X, uniqEV=TRUE) > > > data(iris) > system.time(cN1 <- covNN.1(iris[-5])) user system elapsed 0.09 0.02 0.11 > system.time(cN <- covNNC (iris[-5]))# faster indeed user system elapsed 0.05 0.01 0.06 > > cN2 <- covNNC(iris[-5], trace.lev = 2) nclean.sub(X {150 x 4}, k = 12): alpha.d = 4.9348, k.dists : [1] 0.4003911 0.5316576 0.4729430 0.4358655 0.4630233 0.8165973 0.5080006 [8] 0.4002468 0.7014810 0.4729430 0.5372766 0.4598685 0.4986149 0.6578907 [15] 0.9865180 1.6766405 0.8022084 0.3937506 0.8414153 0.6834719 0.5049454 [22] 0.5015261 0.6893273 0.5153026 0.4875279 0.6019901 0.3574015 0.3824980 [29] 0.4002468 0.4820891 0.4746622 0.4890288 0.9495382 1.2051210 0.4630233 [36] 0.4875279 0.5424999 0.5101600 0.5509906 0.3607570 0.4492088 1.7638065 [43] 0.6024465 0.4978794 0.7328112 0.4986149 0.6834719 0.4986149 0.5015261 [50] 0.4564912 0.8232555 0.6449033 0.6654752 0.8004936 0.5389866 0.5015261 [57] 0.7606229 0.9622246 0.5430063 0.6912283 1.4361149 0.5254180 1.0590853 [64] 0.4596692 0.6214777 0.7147082 0.5509906 0.6101034 1.0553101 0.6840693 [71] 0.7328112 0.5389866 0.7420875 0.5254180 0.5458485 0.5557555 0.7035532 [78] 0.5596444 0.4443544 0.6912283 0.7402335 0.7649075 0.5458485 0.5678788 [85] 0.6834719 0.9746409 0.6078470 0.9594950 0.6298904 0.6905426 0.6180622 [92] 0.5049454 0.5908441 1.0484679 0.5424999 0.5559025 0.4802349 0.5322719 [99] 0.9437765 0.5124970 0.8887088 0.6881182 0.5885882 0.5596444 0.5458485 [106] 0.9880310 1.0198143 0.9413468 0.9843657 1.2753286 0.5804557 0.6064891 [113] 0.5347130 0.8663455 0.7947584 0.5424999 0.5979874 1.8436712 1.5834540 [120] 1.0763259 0.5847012 0.7147082 1.3333053 0.5924970 0.6174851 0.8145461 [127] 0.4961609 0.5131922 0.7004213 0.8026506 1.0190055 1.9358420 0.7328112 [134] 0.5322719 0.7988476 0.9827924 0.9334496 0.5846095 0.5846095 0.5318112 [141] 0.5557555 0.6174851 0.6881182 0.5945271 0.7515754 0.5844867 0.7787172 [148] 0.5318112 0.8527006 0.5783933 EM it=1: loglik.new=5698.06; ll.delta=9.996491e-01, EM it=2: loglik.new=5710.35; ll.delta=2.150854e-03, EM it=3: loglik.new=5717.94; ll.delta=1.328074e-03, EM it=4: loglik.new=5723.28; ll.delta=9.330507e-04, lambda{1,2}= {17.6228,0.953635}; --> z := round(probs): [1] 1 1 1 1 1 0 1 1 1 1 1 1 1 1 0 0 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 [38] 1 1 1 1 0 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1 0 1 0 1 1 1 1 1 0 1 1 1 1 1 [75] 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 0 1 1 1 1 0 1 0 1 1 1 1 0 0 0 0 0 1 [112] 1 1 0 1 1 1 0 0 0 1 1 0 1 1 0 1 1 1 1 0 0 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 [149] 0 1 covNNC(); after nclean*(), --> nnoise= 8 nnoise= 8; ncho1= 4, minv = 4 vv1 vv2 [1,] 0.5284824581307 0.261749007852 [2,] -0.6055151895295 -0.456904892672 [3,] -0.0185844130435 -0.464097859108 [4,] -0.5947371404991 0.712277019924 ot = 1.5707963267948968; extend xcho[,] by 4 rows: k1= 1: pseg= 0.31415926535897937; [sin(.),cos(.)] = [0.30901699437494745, 0.95105651629515353] k1= 2: pseg= 0.62831853071795873; [sin(.),cos(.)] = [0.58778525229247325, 0.80901699437494734] k1= 3: pseg= 0.94247779607693805; [sin(.),cos(.)] = [0.80901699437494745, 0.58778525229247303] k1= 4: pseg= 1.2566370614359175; [sin(.),cos(.)] = [0.95105651629515364, 0.30901699437494728] n2= 8 Mahalanobis D^2: [1] 18.747797 11.163498 17.867061 18.102860 4.349591 5.457991 5.917180 [8] 5.551763 Xa: [1] 13.27670 15.32391 17.37112 19.41833 21.46553 23.51274 -> gap= 2.047208 nclean.sub(X {158 x 4}, k = 12): alpha.d = 4.9348, k.dists : [1] 0.4003911 0.5316576 0.4729430 0.4358655 0.4630233 0.8165973 0.5080006 [8] 0.4002468 0.7014810 0.4729430 0.5372766 0.4598685 0.4986149 0.6578907 [15] 0.9865180 1.6766405 0.8022084 0.3937506 0.8414153 0.6834719 0.5049454 [22] 0.5015261 0.6893273 0.5153026 0.4875279 0.6019901 0.3574015 0.3824980 [29] 0.4002468 0.4820891 0.4746622 0.4890288 0.9495382 1.2051210 0.4630233 [36] 0.4875279 0.5424999 0.5101600 0.5509906 0.3607570 0.4492088 1.6240784 [43] 0.6024465 0.4978794 0.7328112 0.4986149 0.6834719 0.4986149 0.5015261 [50] 0.4564912 0.8232555 0.6449033 0.6654752 0.8004936 0.5389866 0.5015261 [57] 0.7606229 0.9622246 0.5430063 0.6912283 1.4209499 0.5254180 1.0590853 [64] 0.4596692 0.6214777 0.7147082 0.5509906 0.6101034 1.0371046 0.6840693 [71] 0.7328112 0.5389866 0.7420875 0.5254180 0.5458485 0.5557555 0.7035532 [78] 0.5596444 0.4443544 0.6912283 0.7402335 0.7649075 0.5458485 0.5678788 [85] 0.6834719 0.9746409 0.6078470 0.9594950 0.6298904 0.6905426 0.6180622 [92] 0.5049454 0.5908441 1.0484679 0.5424999 0.5559025 0.4802349 0.5322719 [99] 0.9437765 0.5124970 0.8887088 0.6881182 0.5885882 0.5596444 0.5458485 [106] 0.9880310 1.0198143 0.9413468 0.9843657 1.2561456 0.5804557 0.6064891 [113] 0.5347130 0.8663455 0.7947584 0.5424999 0.5979874 1.8300861 1.5834540 [120] 1.0508776 0.5847012 0.7147082 1.3333053 0.5924970 0.6174851 0.8145461 [127] 0.4961609 0.5131922 0.7004213 0.8026506 1.0190055 1.9250383 0.7328112 [134] 0.5322719 0.7988476 0.9827924 0.9334496 0.5846095 0.5846095 0.5318112 [141] 0.5557555 0.6174851 0.6881182 0.5945271 0.7515754 0.5844867 0.7787172 [148] 0.5318112 0.8527006 0.5783933 1.5440475 1.7264589 1.7474427 1.7179325 [155] 1.7807016 1.8284060 1.9362134 2.0435599 EM it=1: loglik.new=5620.36; ll.delta=9.996442e-01, EM it=2: loglik.new=5623.91; ll.delta=6.296881e-04, lambda{1,2}= {16.3654,0.546012}; --> z := round(probs): [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 [38] 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 0 1 1 1 1 1 0 1 1 1 1 1 [75] 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 0 1 1 1 1 0 1 0 1 1 1 1 0 0 0 0 0 1 [112] 1 1 0 1 1 1 0 0 0 1 1 0 1 1 1 1 1 1 1 0 0 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 [149] 0 1 0 0 0 0 0 0 0 0 initv: Named num [1:4] 0.53 0.11 2.823 0.559 - attr(*, "names")= chr [1:4] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width" table( stopv = z := round(probs): [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 [38] 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 0 1 1 1 1 1 0 1 1 1 1 1 [75] 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 0 1 1 1 1 0 1 0 1 1 1 1 0 0 0 0 0 1 [112] 1 1 0 1 1 1 0 0 0 1 1 0 1 1 1 1 1 1 1 0 0 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 [149] 1 1 0 0 0 0 0 0 0 0 time1= 3 table( stopv = z := round(probs): [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 [38] 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 0 1 1 1 1 1 0 1 1 1 1 1 [75] 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 0 1 1 1 1 0 1 0 1 1 1 1 0 0 0 0 0 1 [112] 1 1 1 1 1 1 0 0 0 1 1 0 1 1 1 1 1 1 1 0 0 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 [149] 1 1 0 0 0 0 0 0 0 0 time1= 4 table( stopv = z := round(probs): [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 [38] 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 0 1 1 1 1 1 0 1 1 1 1 1 [75] 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1 0 0 0 0 0 1 [112] 1 1 1 1 1 1 0 0 0 1 1 0 1 1 1 1 1 1 1 0 0 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 [149] 1 1 0 0 0 0 0 0 0 0 time1= 5 table( stopv = z := round(probs): [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 [38] 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 0 1 1 1 1 1 0 1 1 1 1 1 [75] 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1 0 0 0 0 0 1 [112] 1 1 1 1 1 1 0 0 0 1 1 0 1 1 1 1 1 1 1 0 0 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 [149] 1 1 0 0 0 0 0 0 0 0 time1= 6 table( stopv = z := round(probs): [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 [38] 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 0 1 1 1 1 1 0 1 1 1 1 1 [75] 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1 0 0 0 0 0 1 [112] 1 1 1 1 1 1 0 0 0 1 1 0 1 1 1 1 1 1 1 0 0 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 [149] 1 1 0 0 0 0 0 0 0 0 time1= 7 table( stopv = z := round(probs): [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 [38] 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 0 1 1 1 1 1 [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 [112] 1 1 1 1 1 1 0 0 0 1 1 0 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 [149] 1 1 0 0 0 0 0 0 0 0 nclean.sub(X {158 x 4}, k = 12): alpha.d = 4.9348, k.dists : [1] 0.4003911 0.5316576 0.4729430 0.4358655 0.4630233 0.8165973 0.5080006 [8] 0.4002468 0.7014810 0.4729430 0.5372766 0.4598685 0.4986149 0.6578907 [15] 0.9865180 1.6766405 0.8022084 0.3937506 0.8414153 0.6834719 0.5049454 [22] 0.5015261 0.6893273 0.5153026 0.4875279 0.6019901 0.3574015 0.3824980 [29] 0.4002468 0.4820891 0.4746622 0.4890288 0.9495382 1.2051210 0.4630233 [36] 0.4875279 0.5424999 0.5101600 0.5509906 0.3607570 0.4492088 1.6240784 [43] 0.6024465 0.4978794 0.7328112 0.4986149 0.6834719 0.4986149 0.5015261 [50] 0.4564912 0.8232555 0.6449033 0.6654752 0.8004936 0.5389866 0.5015261 [57] 0.7606229 0.9622246 0.5430063 0.6912283 1.4361149 0.5254180 1.0590853 [64] 0.4596692 0.6214777 0.7147082 0.5509906 0.6101034 1.0553101 0.6840693 [71] 0.7328112 0.5389866 0.7420875 0.5254180 0.5458485 0.5557555 0.7035532 [78] 0.5596444 0.4443544 0.6912283 0.7402335 0.7649075 0.5458485 0.5678788 [85] 0.6834719 0.9746409 0.6078470 0.9594950 0.6298904 0.6905426 0.6180622 [92] 0.5049454 0.5908441 1.0484679 0.5424999 0.5559025 0.4802349 0.5322719 [99] 0.9437765 0.5124970 0.8887088 0.6881182 0.5885882 0.5596444 0.5458485 [106] 0.9880310 1.0198143 0.9413468 0.9843657 1.2753286 0.5804557 0.6064891 [113] 0.5347130 0.8663455 0.7947584 0.5424999 0.5979874 1.8300861 1.5834540 [120] 1.0763259 0.5847012 0.7147082 1.3333053 0.5924970 0.6174851 0.8145461 [127] 0.4961609 0.5131922 0.7004213 0.8026506 1.0190055 1.9250383 0.7328112 [134] 0.5322719 0.7988476 0.9827924 0.9334496 0.5846095 0.5846095 0.5318112 [141] 0.5557555 0.6174851 0.6881182 0.5945271 0.7515754 0.5844867 0.7787172 [148] 0.5318112 0.8527006 0.5783933 2.6740900 2.8512315 2.8303116 2.7829071 [155] 2.6054220 2.6630819 2.7874469 2.9025256 EM it=1: loglik.new=4808.34; ll.delta=9.995841e-01, EM it=2: loglik.new=4840.74; ll.delta=6.692834e-03, EM it=3: loglik.new=4872.18; ll.delta=6.451818e-03, EM it=4: loglik.new=4880.52; ll.delta=1.706923e-03, EM it=5: loglik.new=4883.83; ll.delta=6.787204e-04, lambda{1,2}= {10.8317,0.0974727}; --> z := round(probs): [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 [38] 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 0 1 1 1 1 1 [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 [112] 1 1 1 1 1 1 0 0 0 1 1 0 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 [149] 1 1 0 0 0 0 0 0 0 0 > all.equal(cN, cN2, tolerance = 0)# TRUE ? [1] TRUE > stopifnot(all.equal(cN, cN2)) > ## stopifnot(identical(cN, cN2)) > > s1 <- summ.NN(cN1) > ss <- summ.NN(cN) > if(isTRUE(all.equal(ss, s1))) ss else cbind(ss, s1) class pprob incc.p [1,] 1 1.000 1.000 [2,] 1 1.000 1.000 [3,] 1 1.000 1.000 [4,] 1 1.000 1.000 [5,] 1 1.000 1.000 [6,] 1 1.000 0.171 [7,] 1 1.000 1.000 [8,] 1 1.000 1.000 [9,] 1 1.000 1.000 [10,] 1 1.000 1.000 [11,] 1 1.000 1.000 [12,] 1 1.000 1.000 [13,] 1 1.000 1.000 [14,] 1 1.000 1.000 [15,] 1 0.998 0.000 [16,] 0 0.000 0.000 [17,] 1 1.000 0.718 [18,] 1 1.000 1.000 [19,] 1 1.000 0.002 [20,] 1 1.000 1.000 [21,] 1 1.000 1.000 [22,] 1 1.000 1.000 [23,] 1 1.000 1.000 [24,] 1 1.000 1.000 [25,] 1 1.000 1.000 [26,] 1 1.000 1.000 [27,] 1 1.000 1.000 [28,] 1 1.000 1.000 [29,] 1 1.000 1.000 [30,] 1 1.000 1.000 [31,] 1 1.000 1.000 [32,] 1 1.000 1.000 [33,] 1 1.000 0.000 [34,] 0 0.000 0.000 [35,] 1 1.000 1.000 [36,] 1 1.000 1.000 [37,] 1 1.000 1.000 [38,] 1 1.000 1.000 [39,] 1 1.000 1.000 [40,] 1 1.000 1.000 [41,] 1 1.000 1.000 [42,] 0 0.000 0.000 [43,] 1 1.000 1.000 [44,] 1 1.000 1.000 [45,] 1 1.000 1.000 [46,] 1 1.000 1.000 [47,] 1 1.000 1.000 [48,] 1 1.000 1.000 [49,] 1 1.000 1.000 [50,] 1 1.000 1.000 [51,] 1 1.000 0.058 [52,] 1 1.000 1.000 [53,] 1 1.000 1.000 [54,] 1 1.000 0.773 [55,] 1 1.000 1.000 [56,] 1 1.000 1.000 [57,] 1 1.000 0.999 [58,] 1 1.000 0.000 [59,] 1 1.000 1.000 [60,] 1 1.000 1.000 [61,] 0 0.000 0.000 [62,] 1 1.000 1.000 [63,] 0 0.000 0.000 [64,] 1 1.000 1.000 [65,] 1 1.000 1.000 [66,] 1 1.000 1.000 [67,] 1 1.000 1.000 [68,] 1 1.000 1.000 [69,] 0 0.000 0.000 [70,] 1 1.000 1.000 [71,] 1 1.000 1.000 [72,] 1 1.000 1.000 [73,] 1 1.000 1.000 [74,] 1 1.000 1.000 [75,] 1 1.000 1.000 [76,] 1 1.000 1.000 [77,] 1 1.000 1.000 [78,] 1 1.000 1.000 [79,] 1 1.000 1.000 [80,] 1 1.000 1.000 [81,] 1 1.000 1.000 [82,] 1 1.000 0.999 [83,] 1 1.000 1.000 [84,] 1 1.000 1.000 [85,] 1 1.000 1.000 [86,] 1 1.000 0.000 [87,] 1 1.000 1.000 [88,] 1 1.000 0.000 [89,] 1 1.000 1.000 [90,] 1 1.000 1.000 [91,] 1 1.000 1.000 [92,] 1 1.000 1.000 [93,] 1 1.000 1.000 [94,] 0 0.001 0.000 [95,] 1 1.000 1.000 [96,] 1 1.000 1.000 [97,] 1 1.000 1.000 [98,] 1 1.000 1.000 [99,] 1 1.000 0.000 [100,] 1 1.000 1.000 [101,] 1 1.000 0.000 [102,] 1 1.000 1.000 [103,] 1 1.000 1.000 [104,] 1 1.000 1.000 [105,] 1 1.000 1.000 [106,] 1 0.998 0.000 [107,] 0 0.317 0.000 [108,] 1 1.000 0.000 [109,] 1 0.999 0.000 [110,] 0 0.000 0.000 [111,] 1 1.000 1.000 [112,] 1 1.000 1.000 [113,] 1 1.000 1.000 [114,] 1 1.000 0.000 [115,] 1 1.000 0.899 [116,] 1 1.000 1.000 [117,] 1 1.000 1.000 [118,] 0 0.000 0.000 [119,] 0 0.000 0.000 [120,] 0 0.000 0.000 [121,] 1 1.000 1.000 [122,] 1 1.000 1.000 [123,] 0 0.000 0.000 [124,] 1 1.000 1.000 [125,] 1 1.000 1.000 [126,] 1 1.000 0.229 [127,] 1 1.000 1.000 [128,] 1 1.000 1.000 [129,] 1 1.000 1.000 [130,] 1 1.000 0.702 [131,] 0 0.357 0.000 [132,] 0 0.000 0.000 [133,] 1 1.000 1.000 [134,] 1 1.000 1.000 [135,] 1 1.000 0.818 [136,] 1 0.999 0.000 [137,] 1 1.000 0.000 [138,] 1 1.000 1.000 [139,] 1 1.000 1.000 [140,] 1 1.000 1.000 [141,] 1 1.000 1.000 [142,] 1 1.000 1.000 [143,] 1 1.000 1.000 [144,] 1 1.000 1.000 [145,] 1 1.000 1.000 [146,] 1 1.000 1.000 [147,] 1 1.000 0.991 [148,] 1 1.000 1.000 [149,] 1 1.000 0.000 [150,] 1 1.000 1.000 > > try( # testing (tol=0 too small) + chk.NN.new.old(cN, cN1, tol=0) + ) classification accordance matrix: old new 0 1 0 15 0 1 0 135 Error : Not all.equal(UN(cNew[1:4]), UN(cNold[1:4]), tol = tol.1): Component "cov": Mean relative difference: 8.750337e-16 Component "mu": Mean relative difference: 5.177367e-16 Component "postprob": Mean relative difference: 2.279175e-14 > ## This used to fail when we use R's instead of BLAS matrix products: > ## 2026-01: now also "reliably fails" on Linux aa > if(doCheck) + chk.NN.new.old(cN, cN1, tol = 4e-15) # seen 1.1e-15 work classification accordance matrix: old new 0 1 0 15 0 1 0 135 [1] TRUE > > > ## for n = 500, you *do* see it > n <- 500 > set.seed(12) > X <- rbwheel(n, 7, spherize=TRUE) > > lattice::splom(X, cex=.1) > system.time(cNX1 <- covNN.1(X))# 0.82 0.273 user system elapsed 0.42 0.06 0.48 > system.time(cNX <- covNNC (X))# 0.66 0.163 user system elapsed 0.28 0.00 0.28 > system.time(cM <- covMcd (X))# 0.151 0.097 <- ! user system elapsed 0.06 0.00 0.06 > # NB: *slower* times above, when using R's instead of BLAS matrix prod > > try( # --> show relative difference(s): + chk.NN.new.old(cNX, cNX1, tol=0) + ) classification accordance matrix: old new 0 1 0 179 0 1 0 321 Error : Not all.equal(UN(cNew[1:4]), UN(cNold[1:4]), tol = tol.1): Component "cov": Mean relative difference: 3.901356e-15 Component "mu": Mean relative difference: 1.538281e-14 Component "postprob": Mean relative difference: 7.168447e-15 > if(doCheck && ourBLAS) # did fail with ATLAS in R-devel 2023-1-1 + chk.NN.new.old(cNX, cNX1) > > stopifnot(exprs = { + all.equal(1900.4208, print(kappa(cM $cov)))# 1990.8.. then 1900.421 + all.equal(4.437255 , print(kappa(cNX$cov)), tolerance = 5e-3)# was 4.44858 + all.equal(1.04778125, print(kappa(cov(X)))) + }) [1] 1900.421 [1] 4.437255 [1] 1.047781 > > ## ---- d = 1 : > X1 <- cbind(c(1:6, 1000)) > > var(X1) [,1] [1,] 141861.8 > ## [,1] > ## [1,] 141861.8 > ## if 1000 was not an outlier: > var(1:7) ## 4.666667 [1] 4.666667 > > covNNC(X1)$cov ## -- really not at all robust: [,1] [1,] 121595.8 > ## [,1] > ## [1,] 121595.8 > > (C.mcd <- covMcd(X1)$cov) [,1] [1,] 4.884832 > Cm <- as.matrix(if(newRB) 4.8848 else 7.790004) > all.equal(Cm, C.mcd, tolerance=0) # 6.633e-6 [1] "Mean relative difference: 6.633081e-06" > stopifnot(all.equal(Cm, C.mcd, tol = 2e-5)) > > > MASS::cov.rob(X1)$cov [,1] [1,] 3.5 > ## [,1] > ## [1,] 3.5 > (C.B <- BACON(X1)$cov) rank(x.ord[1:m,] >= p ==> chosen m = 3 MV-BACON (subset no. 1): 3 of 7 (42.86 %) MV-BACON (subset no. 2): 6 of 7 (85.71 %) MV-BACON (subset no. 3): 6 of 7 (85.71 %) [,1] [1,] 3.5 > ## [,1] > ## [1,] 3.5 > all.equal(C.B, as.matrix(3.5), tol=0) [1] TRUE > stopifnot(all.equal(C.B, as.matrix(3.5))) > > if(FALSE) ## FIXME (in robustbase!): should work for p=1 + covOGK(X1)$cov > > > ## Less trivial data --- also used in ../man/BACON.Rd > data(starsCYG, package = "robustbase") > > op <- options(warn = 2)# no warnings allowed > str(B.ST <- with(starsCYG, BACON(x = log.Te, y = log.light))) rank(x.ord[1:m,] >= p ==> chosen m = 4 MV-BACON (subset no. 1): 4 of 47 (8.51 %) MV-BACON (subset no. 2): 5 of 47 (10.64 %) MV-BACON (subset no. 3): 5 of 47 (10.64 %) Reg-BACON (init subset no. 0): 8 of 47 (17.02 %) Reg-BACON (init subset no. 0): 3 of 47 (6.38 %) Reg-BACON (init subset no. 1): 4 of 47 (8.51 %) Reg-BACON (init subset no. 2): 5 of 47 (10.64 %) Reg-BACON (init subset no. 3): 6 of 47 (12.77 %) Reg-BACON (init subset no. 4): 7 of 47 (14.89 %) Reg-BACON (init subset no. 5): 8 of 47 (17.02 %) Reg-BACON (subset no. 1): 15 of 47 (31.91 %) Reg-BACON (subset no. 2): 28 of 47 (59.57 %) Reg-BACON (subset no. 3): 40 of 47 (85.11 %) Reg-BACON (subset no. 4): 42 of 47 (89.36 %) Reg-BACON (subset no. 5): 43 of 47 (91.49 %) Reg-BACON (subset no. 6): 43 of 47 (91.49 %) List of 5 $ subset : logi [1:47] TRUE TRUE TRUE TRUE TRUE TRUE ... $ tis : num [1:47] 0.854 1.176 0.674 1.176 1.116 ... $ mv.subset: logi [1:47] FALSE FALSE FALSE FALSE TRUE FALSE ... $ mv.dis : num [1:47] 17.44 59.93 7.16 59.93 1.79 ... $ steps : Named int [1:2] 3 6 ..- attr(*, "names")= chr [1:2] "mv" "lm" > Bgood <- which(B.ST$subset) > cat("starsCYG - outliers (BACON):"); setdiff(seq_len(nrow(starsCYG)), Bgood) starsCYG - outliers (BACON):[1] 11 20 30 34 > .Platform$r_arch # typical empty for Linux [1] "x64" > ## 32-bit <-> 64-bit different results {tested on Linux & Windows Server} > is32 <- .Machine$sizeof.pointer == 4 ## <- should work for Linux/MacOS/Windows > isMac <- Sys.info()[["sysname"]] == "Darwin" > isSun <- Sys.info()[["sysname"]] == "SunOS" > isWin <- .Platform$OS.type == "windows" > (Platf_arch <- paste0(.Platform$r_arch, # maybe distinguish more, e.g. "no-ldouble", "BLAS" version , ?? + if(isWin) "_win")) [1] "x64_win" > knownPl <- Platf_arch %in% c("x64", "i386_win", "i386") # update! > stopifnot(exprs = { + switch( + Platf_arch, + "i386_win" =, # Platform: i386-w64-mingw32/i386 (32-bit); Windows Server x64 (build 14393) + "x64" = identical(Bgood, c(25L, 27:29, 33L, 38L, 43L, 45L)), + "i386" = identical(Bgood, c(25L, 27:29, 33L, 35:36, 38L, 42L, 43L, 45L)), + # older version of i386-windows + ## other platforms: + { + message("Platform architecture (see above) not yet tested for BACON result") + TRUE + }) + }) Platform architecture (see above) not yet tested for BACON result > > > plot(starsCYG) > lmST <- lm(log.light ~ log.Te, data = starsCYG) > abline(lmST, col = "gray") # least squares line > ## 'subset': A good set of of points (to determine regression): > colB <- adjustcolor(2, 1/2) > points(log.light ~ log.Te, data = starsCYG, subset = B.ST$subset, + pch = 19, cex = 1.5, col = colB) > ## A BACON-derived line: > lmB <- lm(log.light ~ log.Te, data = starsCYG, subset = B.ST$subset) > abline(lmB, col = colB, lwd = 2) > cfT <- switch(Platf_arch # use lm(..., subset = sfsmisc::inverseWhich(Bgood)) + , "x64" = c(-10.130429, 3.4450151) + , "i386" = c(-9.5759001, 3.3109007) + ## otherwise: + , NULL) > cf <- unname(coef(lmB)) > dput(signif(cf, 8)) c(-4.0565237, 2.0466574) > if(!is.null(cfT)) withAutoprint({ + all.equal(cf, cfT, tolerance = 0)# 64b: 6.4e-10 + stopifnot(all.equal(cf, cfT)) + }) > > options(op) > > proc.time() user system elapsed 2.18 0.21 2.39