R Under development (unstable) (2024-08-16 r87026 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 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.

> dodata <- function(nrep=1, time=FALSE, short=FALSE, full=TRUE, method = c("FASTMVE","MASS")){
+ ##@bdescr
+ ## Test the function covMve() on the literature datasets:
+ ##
+ ## Call covMve() for all regression datasets available in rrco/robustbasev and print:
+ ##  - execution time (if time == TRUE)
+ ##  - objective fucntion
+ ##  - best subsample found (if short == false)
+ ##  - outliers identified (with cutoff 0.975) (if short == false)
+ ##  - estimated center and covarinance matrix if full == TRUE)
+ ##
+ ##@edescr
+ ##
+ ##@in  nrep              : [integer] number of repetitions to use for estimating the
+ ##                                   (average) execution time
+ ##@in  time              : [boolean] whether to evaluate the execution time
+ ##@in  short             : [boolean] whether to do short output (i.e. only the
+ ##                                   objective function value). If short == FALSE,
+ ##                                   the best subsample and the identified outliers are
+ ##                                   printed. See also the parameter full below
+ ##@in  full              : [boolean] whether to print the estimated cente and covariance matrix
+ ##@in  method            : [character] select a method: one of (FASTMCD, MASS)
+ 
+     domve <- function(x, xname, nrep=1){
+         n <- dim(x)[1]
+         p <- dim(x)[2]
+         alpha <- 0.5
+         h <- h.alpha.n(alpha, n, p)
+         if(method == "MASS"){
+             mve <- cov.mve(x, quantile.used=h)
+             quan <- h   #default: floor((n+p+1)/2)
+             crit <- mve$crit
+             best <- mve$best
+             mah <- mahalanobis(x, mve$center, mve$cov)
+             quantiel <- qchisq(0.975, p)
+             wt <- as.numeric(mah < quantiel)
+         }
+         else{
+             mve <- CovMve(x, trace=FALSE)
+             quan <- as.integer(mve@quan)
+             crit <- log(mve@crit)
+             best <- mve@best
+             wt <- mve@wt
+         }
+ 
+ 
+         if(time){
+            xtime <- system.time(dorep(x, nrep, method))[1]/nrep
+            xres <- sprintf("%3d %3d %3d %12.6f %10.3f\n", dim(x)[1], dim(x)[2], quan, crit, xtime)
+         }
+         else{
+             xres <- sprintf("%3d %3d %3d %12.6f\n", dim(x)[1], dim(x)[2], quan, crit)
+         }
+ 
+         lpad<-lname-nchar(xname)
+         cat(pad.right(xname,lpad), xres)
+ 
+         if(!short){
+             cat("Best subsample: \n")
+             print(best)
+ 
+             ibad <- which(wt == 0)
+             names(ibad) <- NULL
+             nbad <- length(ibad)
+             cat("Outliers: ", nbad, "\n")
+             if(nbad > 0)
+                 print(ibad)
+             if(full){
+                 cat("-------------\n")
+                 show(mve)
+             }
+             cat("--------------------------------------------------------\n")
+         }
+     }
+ 
+     options(digits = 5)
+     set.seed(101) # <<-- sub-sampling algorithm now based on R's RNG and seed
+ 
+     lname <- 20
+ 
+     ## VT::15.09.2013 - this will render the output independent
+     ##  from the version of the package
+     suppressPackageStartupMessages(library(rrcov))
+ 
+     method <- match.arg(method)
+     if(method == "MASS")
+         library(MASS)
+ 
+ 
+     data(heart)
+     data(starsCYG)
+     data(phosphor)
+     data(stackloss)
+     data(coleman)
+     data(salinity)
+     data(wood)
+ 
+     data(hbk)
+ 
+     data(Animals, package = "MASS")
+     brain <- Animals[c(1:24, 26:25, 27:28),]
+     data(milk)
+     data(bushfire)
+ 
+     tmp <- sys.call()
+     cat("\nCall: ", deparse(substitute(tmp)),"\n")
+ 
+     cat("Data Set               n   p  Half LOG(obj)        Time\n")
+     cat("========================================================\n")
+     domve(heart[, 1:2], data(heart), nrep)
+     domve(starsCYG, data(starsCYG), nrep)
+     domve(data.matrix(subset(phosphor, select = -plant)), data(phosphor), nrep)
+     domve(stack.x, data(stackloss), nrep)
+     domve(data.matrix(subset(coleman, select = -Y)), data(coleman), nrep)
+     domve(data.matrix(subset(salinity, select = -Y)), data(salinity), nrep)
+     domve(data.matrix(subset(wood, select = -y)), data(wood), nrep)
+     domve(data.matrix(subset(hbk,  select = -Y)),data(hbk), nrep)
+ 
+     domve(brain, "Animals", nrep)
+     domve(milk, data(milk), nrep)
+     domve(bushfire, data(bushfire), nrep)
+     cat("========================================================\n")
+ }
> 
> dogen <- function(nrep=1, eps=0.49, method=c("FASTMVE", "MASS")){
+ 
+     domve <- function(x, nrep=1){
+         gc()
+         xtime <- system.time(dorep(x, nrep, method))[1]/nrep
+         cat(sprintf("%6d %3d %10.2f\n", dim(x)[1], dim(x)[2], xtime))
+         xtime
+     }
+ 
+     set.seed(1234)
+ 
+     ## VT::15.09.2013 - this will render the output independent
+     ##  from the version of the package
+     suppressPackageStartupMessages(library(rrcov))
+     library(MASS)
+ 
+     method <- match.arg(method)
+ 
+     ap <- c(2, 5, 10, 20, 30)
+     an <- c(100, 500, 1000, 10000, 50000)
+ 
+     tottime <- 0
+     cat("     n   p       Time\n")
+     cat("=====================\n")
+     for(i in 1:length(an)) {
+         for(j in 1:length(ap)) {
+             n <- an[i]
+             p <- ap[j]
+             if(5*p <= n){
+                 xx <- gendata(n, p, eps)
+                 X <- xx$X
+                 tottime <- tottime + domve(X, nrep)
+             }
+         }
+     }
+ 
+     cat("=====================\n")
+     cat("Total time: ", tottime*nrep, "\n")
+ }
> 
> docheck <- function(n, p, eps){
+     xx <- gendata(n,p,eps)
+     mve <- CovMve(xx$X)
+     check(mve, xx$xind)
+ }
> 
> check <- function(mcd, xind){
+ ##  check if mcd is robust w.r.t xind, i.e. check how many of xind
+ ##  did not get zero weight
+     mymatch <- xind %in% which(mcd@wt == 0)
+     length(xind) - length(which(mymatch))
+ }
> 
> dorep <- function(x, nrep=1, method=c("FASTMVE","MASS")){
+ 
+     method <- match.arg(method)
+     for(i in 1:nrep)
+     if(method == "MASS")
+         cov.mve(x)
+     else
+         CovMve(x)
+ }
> 
> #### gendata() ####
> # Generates a location contaminated multivariate
> # normal sample of n observations in p dimensions
> #    (1-eps)*Np(0,Ip) + eps*Np(m,Ip)
> # where
> #    m = (b,b,...,b)
> # Defaults: eps=0 and b=10
> #
> gendata <- function(n,p,eps=0,b=10){
+ 
+     if(missing(n) || missing(p))
+         stop("Please specify (n,p)")
+     if(eps < 0 || eps >= 0.5)
+         stop(message="eps must be in [0,0.5)")
+     X <- mvrnorm(n,rep(0,p),diag(1,nrow=p,ncol=p))
+     nbad <- as.integer(eps * n)
+     if(nbad > 0){
+         Xbad <- mvrnorm(nbad,rep(b,p),diag(1,nrow=p,ncol=p))
+         xind <- sample(n,nbad)
+         X[xind,] <- Xbad
+     }
+     list(X=X, xind=xind)
+ }
> 
> pad.right <- function(z, pads)
+ {
+ ### Pads spaces to right of text
+     padding <- paste(rep(" ", pads), collapse = "")
+     paste(z, padding, sep = "")
+ }
> 
> whatis<-function(x){
+     if(is.data.frame(x))
+         cat("Type: data.frame\n")
+     else if(is.matrix(x))
+         cat("Type: matrix\n")
+     else if(is.vector(x))
+         cat("Type: vector\n")
+     else
+         cat("Type: don't know\n")
+ }
> 
> ## VT::15.09.2013 - this will render the output independent
> ##  from the version of the package
> suppressPackageStartupMessages(library(rrcov))
> 
> dodata()

Call:  dodata() 
Data Set               n   p  Half LOG(obj)        Time
========================================================
heart                 12   2   7     3.827606
Best subsample: 
[1]  1  4  7  8  9 10 11
Outliers:  3 
[1]  2  6 12
-------------

Call:
CovMve(x = x, trace = FALSE)
-> Method:  Minimum volume ellipsoid estimator 

Robust Estimate of Location: 
height  weight  
  34.9    27.0  

Robust Estimate of Covariance: 
        height  weight
height  142     217   
weight  217     350   
--------------------------------------------------------
starsCYG              47   2  25    -2.742997
Best subsample: 
 [1]  2  4  6  8 12 13 16 23 24 25 26 28 31 32 33 37 38 39 41 42 43 44 45 46 47
Outliers:  7 
[1]  7  9 11 14 20 30 34
-------------

Call:
CovMve(x = x, trace = FALSE)
-> Method:  Minimum volume ellipsoid estimator 

Robust Estimate of Location: 
   log.Te  log.light  
     4.41       4.93  

Robust Estimate of Covariance: 
           log.Te  log.light
log.Te     0.0173  0.0578   
log.light  0.0578  0.3615   
--------------------------------------------------------
phosphor              18   2  10     4.443101
Best subsample: 
 [1]  3  5  8  9 11 12 13 14 15 17
Outliers:  3 
[1]  1  6 10
-------------

Call:
CovMve(x = x, trace = FALSE)
-> Method:  Minimum volume ellipsoid estimator 

Robust Estimate of Location: 
  inorg  organic  
   15.2     39.4  

Robust Estimate of Covariance: 
         inorg  organic
inorg    188    230    
organic  230    339    
--------------------------------------------------------
stackloss             21   3  12     3.327582
Best subsample: 
 [1]  4  5  6  7  8  9 10 11 12 13 14 20
Outliers:  3 
[1] 1 2 3
-------------

Call:
CovMve(x = x, trace = FALSE)
-> Method:  Minimum volume ellipsoid estimator 

Robust Estimate of Location: 
  Air.Flow  Water.Temp  Acid.Conc.  
      56.7        20.2        85.5  

Robust Estimate of Covariance: 
            Air.Flow  Water.Temp  Acid.Conc.
Air.Flow    34.31     11.07       23.54     
Water.Temp  11.07      9.23        7.85     
Acid.Conc.  23.54      7.85       47.35     
--------------------------------------------------------
coleman               20   5  13     2.065143
Best subsample: 
 [1]  1  3  4  5  7  8 11 14 16 17 18 19 20
Outliers:  5 
[1]  2  6  9 10 13
-------------

Call:
CovMve(x = x, trace = FALSE)
-> Method:  Minimum volume ellipsoid estimator 

Robust Estimate of Location: 
  salaryP   fatherWc    sstatus  teacherSc  motherLev  
     2.79      44.26       3.59      25.08       6.38  

Robust Estimate of Covariance: 
           salaryP   fatherWc  sstatus   teacherSc  motherLev
salaryP      0.2920    1.1188    2.0421    0.3487     0.0748 
fatherWc     1.1188  996.7540  338.6587    7.1673    23.1783 
sstatus      2.0421  338.6587  148.2501    4.4894     7.8135 
teacherSc    0.3487    7.1673    4.4894    0.9082     0.3204 
motherLev    0.0748   23.1783    7.8135    0.3204     0.6024 
--------------------------------------------------------
salinity              28   3  16     2.002555
Best subsample: 
 [1]  1  7  8  9 12 13 14 18 19 20 21 22 25 26 27 28
Outliers:  5 
[1]  5 11 16 23 24
-------------

Call:
CovMve(x = x, trace = FALSE)
-> Method:  Minimum volume ellipsoid estimator 

Robust Estimate of Location: 
  X1    X2    X3  
10.2   3.1  22.4  

Robust Estimate of Covariance: 
    X1      X2      X3    
X1  14.387   1.153  -4.072
X2   1.153   5.005  -0.954
X3  -4.072  -0.954   2.222
--------------------------------------------------------
wood                  20   5  13    -5.471407
Best subsample: 
 [1]  1  2  3  5  9 10 12 13 14 15 17 18 20
Outliers:  5 
[1]  4  6  8 11 19
-------------

Call:
CovMve(x = x, trace = FALSE)
-> Method:  Minimum volume ellipsoid estimator 

Robust Estimate of Location: 
   x1     x2     x3     x4     x5  
0.576  0.123  0.531  0.538  0.889  

Robust Estimate of Covariance: 
    x1         x2         x3         x4         x5       
x1   7.45e-03   1.11e-03   1.83e-03  -2.90e-05  -5.65e-04
x2   1.11e-03   3.11e-04   7.68e-04   3.37e-05   3.85e-05
x3   1.83e-03   7.68e-04   4.30e-03  -9.96e-04  -6.27e-05
x4  -2.90e-05   3.37e-05  -9.96e-04   3.02e-03   1.91e-03
x5  -5.65e-04   3.85e-05  -6.27e-05   1.91e-03   2.25e-03
--------------------------------------------------------
hbk                   75   3  39     1.096831
Best subsample: 
 [1] 15 17 18 19 20 21 24 27 28 30 32 33 35 36 40 41 42 43 44 46 48 49 50 53 54
[26] 55 56 58 59 64 65 66 67 70 71 72 73 74 75
Outliers:  14 
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14
-------------

Call:
CovMve(x = x, trace = FALSE)
-> Method:  Minimum volume ellipsoid estimator 

Robust Estimate of Location: 
  X1    X2    X3  
1.48  1.86  1.73  

Robust Estimate of Covariance: 
    X1     X2     X3   
X1  1.695  0.230  0.265
X2  0.230  1.679  0.119
X3  0.265  0.119  1.683
--------------------------------------------------------
Animals               28   2  15     8.945423
Best subsample: 
 [1]  1  3  4  5 10 11 17 18 21 22 23 24 26 27 28
Outliers:  9 
[1]  2  6  7  9 12 14 15 16 25
-------------

Call:
CovMve(x = x, trace = FALSE)
-> Method:  Minimum volume ellipsoid estimator 

Robust Estimate of Location: 
 body  brain  
 48.3  127.3  

Robust Estimate of Covariance: 
       body   brain
body   10767  16872
brain  16872  46918
--------------------------------------------------------
milk                  86   8  47    -1.160085
Best subsample: 
 [1]  4  5  7  8  9 10 11 19 21 22 23 24 26 30 31 33 34 35 36 37 38 39 42 43 45
[26] 46 54 56 57 59 60 61 62 63 64 65 66 67 69 72 76 78 79 81 82 83 85
Outliers:  18 
 [1]  1  2  3 12 13 14 15 16 17 18 20 27 41 44 47 70 74 75
-------------

Call:
CovMve(x = x, trace = FALSE)
-> Method:  Minimum volume ellipsoid estimator 

Robust Estimate of Location: 
    X1      X2      X3      X4      X5      X6      X7      X8  
  1.03   35.91   33.02   26.08   25.06   24.99  122.93   14.38  

Robust Estimate of Covariance: 
    X1        X2        X3        X4        X5        X6        X7      
X1  6.00e-07  1.51e-04  3.34e-04  3.09e-04  2.82e-04  2.77e-04  1.09e-03
X2  1.51e-04  2.03e+00  3.83e-01  3.04e-01  2.20e-01  3.51e-01  2.18e+00
X3  3.34e-04  3.83e-01  1.58e+00  1.21e+00  1.18e+00  1.20e+00  1.60e+00
X4  3.09e-04  3.04e-01  1.21e+00  9.82e-01  9.39e-01  9.53e-01  1.36e+00
X5  2.82e-04  2.20e-01  1.18e+00  9.39e-01  9.67e-01  9.52e-01  1.34e+00
X6  2.77e-04  3.51e-01  1.20e+00  9.53e-01  9.52e-01  9.92e-01  1.38e+00
X7  1.09e-03  2.18e+00  1.60e+00  1.36e+00  1.34e+00  1.38e+00  6.73e+00
X8  3.33e-05  2.92e-01  2.65e-01  1.83e-01  1.65e-01  1.76e-01  5.64e-01
    X8      
X1  3.33e-05
X2  2.92e-01
X3  2.65e-01
X4  1.83e-01
X5  1.65e-01
X6  1.76e-01
X7  5.64e-01
X8  1.80e-01
--------------------------------------------------------
bushfire              38   5  22     5.644315
Best subsample: 
 [1]  1  2  3  4  5  6 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
Outliers:  15 
 [1]  7  8  9 10 11 29 30 31 32 33 34 35 36 37 38
-------------

Call:
CovMve(x = x, trace = FALSE)
-> Method:  Minimum volume ellipsoid estimator 

Robust Estimate of Location: 
 V1   V2   V3   V4   V5  
107  147  263  215  277  

Robust Estimate of Covariance: 
    V1     V2     V3     V4     V5   
V1    519    375  -2799   -619   -509
V2    375    320  -1671   -342   -289
V3  -2799  -1671  18373   4314   3480
V4   -619   -342   4314   1076    854
V5   -509   -289   3480    854    683
--------------------------------------------------------
========================================================
> 
> proc.time()
   user  system elapsed 
   0.59    0.12    0.70