R Under development (unstable) (2024-08-15 r87022 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.

> #### Tests psi(), chi(),... etc and  tuning.psi, tuning.chi :
> 
> library(robustbase)
> source(system.file("xtraR/plot-psiFun.R", package = "robustbase", mustWork=TRUE))
> source(system.file("xtraR/test-tools.R",  package = "robustbase")) # assert.EQ
> ## NB: Code to compute some of the constants/tuning parameters (e.g. in ../R/lmrob.MM.R ):
> ## --  MM's ~/R/MM/Pkg-ex/robustbase/lmrob-check-const.R and .../psi-tuning-hampel.R
> 
> (doExtras <- robustbase:::doExtras())
[1] FALSE
> 
> ### (1) Test the functions themselves --------------------------------
> if(!dev.interactive(orNone=TRUE)) pdf("rob-psifns.pdf")
> 
> ## Simple version, no error checking, no derivative, nothing:
> psiGGW <- function(x, a,b,c) {
+     ifelse((ax <- abs(x)) < c,
+            x,
+            ifelse((ea <- -((ax-c)^b)/(2*a)) < -708.4, 0, x * exp(ea)))
+ }
> assert.EQ(Mpsi  (5:9, cc=c(0, a=1/8,b=2,c=1/8, NA), "GGW"),
+           psiGGW(5:9,	      a=1/8,b=2,c=1/8), tol = 1e-13)
> 
> 
> ## Check that psi(<empty>)  |->  <empty>  works; ditto for +-Inf, NA,..
> cG <- c(-.5, 1, .95, NA) # one of the 6 "builtin"s
> d0 <- numeric()
> IoI <- c(-Inf, 0, Inf)
> NN <- c(NaN, NA)
> 
> cGs <- list(  c(-.4, 1.5,    0.85,  NA)
+             , c(-.4, 1.5 ,   0.90,  NA)
+             , c(-.4, 1.5 ,   0.95,  NA)
+             , c(-.4, 1.5,    0.975, NA)
+             , c(-.4, 1.5,    0.99 , NA)
+             , c(-.4, 1.5,    0.995, NA)
+             ##
+             , c(-.4, 1.25,   0.975, NA)
+             , c(-.4, 1.1,    0.975, NA)
+             , c(-.4, 1.025,  0.975, NA)
+             , c(-.4, 1.0125, 0.975, NA)
+             ##
+             ## FIXME , c(-.1, 1.25, 0.95, NA)
+             ## FIXME , c(-.1, 1.25, 0.99, NA)
+             )
> st <- system.time(
+ cG.cnst <- lapply(cGs, function(cc)
+                   lmrob.control(psi = "ggw", tuning.psi = cc)$tuning.psi)
+ )
> cat('Time for constants computation of tuning.psi: ', st,'\n')
Time for constants computation of tuning.psi:  0.27 0 0.27 NA NA 
> cGct <- t(sapply(cG.cnst, attr, "constants"))[,-1]
> colnames(cGct) <- c("a","b","c", "rhoInf")
> signif(cGct, 4)
           a     b      c rhoInf
 [1,] 1.0170 1.500 0.4996  2.384
 [2,] 1.2810 1.500 0.5826  3.242
 [3,] 1.8100 1.500 0.7335  5.139
 [4,] 2.4430 1.500 0.8959  7.666
 [5,] 3.4380 1.500 1.1250 12.090
 [6,] 4.2970 1.500 1.3050 16.280
 [7,] 1.3780 1.250 1.4350  7.654
 [8,] 1.0140 1.100 1.7000  7.643
 [9,] 0.8873 1.025 1.8130  7.712
[10,] 0.8693 1.012 1.8300  7.733
> assert.EQ(sapply(cG.cnst, function(cc) MrhoInf(cc, "ggw")),
+           cGct[,"rhoInf"], tol = 1e-8)
> 
> 
> ## Do these checks for a *list* of (c.par, psi) combinations:
> c.psi.list <- list(
+     list(1.345, "Huber"),
+     list(1.8,   "Huber"),
+     list(cG, "GGW"),
+     list(c(2,4,8), "Hampel"),
+     list(c(1.5,3.5,8)*0.90, "Hampel"),
+     list(par=c(-.5,1.5,.95,NA), "lqq"),
+     list(bcs=c(1, 1, 1.25), "lqq"),
+     list(1.1, "optimal"),
+     list(0.1, "optimal"),
+     list(2.3, "Welsh")
+     )
> 
> for(c.psi in c.psi.list) {
+     tPar <-  c.psi[[1]]; psi <- c.psi[[2]]
+     stopifnot(is.numeric(tPar), is.character(psi))
+     cat("Psi function ", psi,"; tuning par. c[]= (",
+         paste(formatC(tPar, width=1), collapse=", "),")\n")
+     for(FUN in list(Mpsi, Mchi, Mwgt))
+ 	stopifnot(identical(d0, FUN(d0, tPar, psi=psi)),
+                   identical(NN, FUN(NN, tPar, psi=psi)))
+     stopifnot(identical(c(0,1,0), Mwgt(IoI, tPar,psi=psi)))
+     if(isPsi.redesc(psi))
+ 	stopifnot(identical(c(0,0,0), Mpsi(IoI, tPar,psi=psi)),
+ 		  identical(c(1,0,1), Mchi(IoI, tPar,psi=psi)))
+     else if(psi == "Huber") {
+ 	stopifnot(identical(c(-tPar,0,tPar), Mpsi(IoI, tPar,psi=psi)),
+ 		  identical(c(  Inf,0, Inf), Mchi(IoI, tPar,psi=psi)))
+     }
+     cat("chkPsi..(): ")
+     isHH <- psi %in% c("Huber", "Hampel") # not differentiable
+     tol <- switch(tolower(psi),
+                   "huber"=, "hampel"= c(.001, 1.0),
+                   "optimal" = .008,
+                   "ggw" = c(5e-5, 5e-3, 1e-12),
+                   "lqq" = c(1e-5, 5e-5, 1e-5, .08)) # .08 needed for bcs=c(1, 1, 1.25)
+     if(is.null(tol)) tol <- 1e-4 # default otherwise
+     cc <- chkPsi..(c(-5, 10), psi=psi, par=tPar, doD2 = !isHH, tol=tol)
+     ##    --------
+     cc. <- cc[!is.na(cc)]
+     if(is.logical(cc) && all(cc.))
+ 	cat(" [Ok]\n")
+     else {
+ 	cat(" not all Ok:\n")
+ 	print(cc.[cc. != "TRUE"])
+     }
+     cat("------------------------\n\n")
+ }
Psi function  Huber ; tuning par. c[]= ( 1.345 )
chkPsi..(): Not checking psi''() := Mpsi(*, deriv=2)
 [Ok]
------------------------

Psi function  Huber ; tuning par. c[]= ( 1.8 )
chkPsi..(): Not checking psi''() := Mpsi(*, deriv=2)
 [Ok]
------------------------

Psi function  GGW ; tuning par. c[]= ( -0.5, 1, 0.95, NA )
chkPsi..(): Not checking psi''() := Mpsi(*, deriv=2)
 [Ok]
------------------------

Psi function  Hampel ; tuning par. c[]= ( 2, 4, 8 )
chkPsi..(): Not checking psi''() := Mpsi(*, deriv=2)
 [Ok]
------------------------

Psi function  Hampel ; tuning par. c[]= ( 1.35, 3.15, 7.2 )
chkPsi..(): Not checking psi''() := Mpsi(*, deriv=2)
 [Ok]
------------------------

Psi function  lqq ; tuning par. c[]= ( -0.5, 1.5, 0.95, NA )
chkPsi..():  [Ok]
------------------------

Psi function  lqq ; tuning par. c[]= ( 1, 1, 1.25 )
chkPsi..():  [Ok]
------------------------

Psi function  optimal ; tuning par. c[]= ( 1.1 )
chkPsi..(): Not checking psi''() := Mpsi(*, deriv=2)
 [Ok]
------------------------

Psi function  optimal ; tuning par. c[]= ( 0.1 )
chkPsi..(): Not checking psi''() := Mpsi(*, deriv=2)
 [Ok]
------------------------

Psi function  Welsh ; tuning par. c[]= ( 2.3 )
chkPsi..(): Not checking psi''() := Mpsi(*, deriv=2)
 [Ok]
------------------------

> 
> ## Demonstrate lmrob.tau.fast.coefs() computations (well, to 5--7 digits accuracy):
> ## a) values from the if(fast && ..)  switch(control$psi, .....) in lmrob.tau() { --> ../R/lmrob.MM.R }
> lmrob.tauFC <- list(optimal = list(c = 1.060158, tfc = c(tfact = 0.94735878, tcorr = -0.09444537)),
+                     bisquare= list(c = 4.685061, tfc = c(tfact = 0.9473684 , tcorr = -0.0900833 )),
+                     welsh   = list(c = 2.11,     tfc = c(tfact = 0.94732953, tcorr = -0.07569506)),
+                     ggw_1   = list(c = c(-.5, 1.0, 0.95, NA),     tfc = c(tfact = 0.9473787 , tcorr = -0.1143846)),
+                     ggw_2   = list(c = c(-.5, 1.5, 0.95, NA),     tfc = c(tfact = 0.94741036, tcorr = -0.08424648)),
+                     lqq     = list(c = c(-.5, 1.5, 0.95, NA),     tfc = c(tfact = 0.94736359, tcorr = -0.08594805)),
+                     hampel  = list(c = c(1.5, 3.5, 8)* 0.9016085, tfc = c(tfact = 0.94739770, tcorr = -0.04103958)))
> lmrob.tau.fast.coefs <- robustbase:::lmrob.tau.fast.coefs
> rr <- sapply(names(lmrob.tauFC), function(nm) { L <- lmrob.tauFC[[nm]]
+     tfc <- lmrob.tau.fast.coefs(cc = L$c, psi = sub("_[0-9]*$", '', nm))
+     ## workaround (for older robustbase:) names(tfc)[2] <- "tcorr"
+     Dif <- all.equal(L$tfc, tfc, tolerance = 0)
+     cat(nm,": 'true' vs lmrob.tau.fast.coefs()-computed: ", Dif, "\n")
+     stopifnot(all.equal(L$tfc, tfc, tol = 4e-5)) # default tol 1.49e-8  is way too small
+     if(doExtras) { ## compute "slow coefs"
+         ## Signature:  lmrob.tau(obj, x=obj$x, control = obj$control, h, fast=TRUE)
+         ##   h = hatvalues  {are made from 'obj' if missing}
+ ### TODO
+         if(FALSE) tau <- lmrob.tau(list(), x=, control=lmrob.control(psi = nm ), fast=FALSE)
+         ## .. and now ?
+     }
+     c(tfc, rel.diff = as.numeric(sub(".*:", '', Dif)))
+     })
optimal : 'true' vs lmrob.tau.fast.coefs()-computed:  Mean relative difference: 7.030091e-07 
bisquare : 'true' vs lmrob.tau.fast.coefs()-computed:  Mean relative difference: 1.127089e-07 
welsh : 'true' vs lmrob.tau.fast.coefs()-computed:  Mean relative difference: 7.957021e-09 
ggw_1 : 'true' vs lmrob.tau.fast.coefs()-computed:  Mean relative difference: 5.127632e-06 
ggw_2 : 'true' vs lmrob.tau.fast.coefs()-computed:  Mean relative difference: 1.814587e-05 
lqq : 'true' vs lmrob.tau.fast.coefs()-computed:  Mean relative difference: 1.722337e-05 
hampel : 'true' vs lmrob.tau.fast.coefs()-computed:  Mean relative difference: 3.291664e-05 
> t(rr)
             tfact       tcorr     rel.diff
optimal  0.9473584 -0.09444571 7.030091e-07
bisquare 0.9473684 -0.09008338 1.127089e-07
welsh    0.9473295 -0.07569506 7.957021e-09
ggw_1    0.9473812 -0.11438166 5.127632e-06
ggw_2    0.9474191 -0.08423651 1.814587e-05
lqq      0.9473721 -0.08593877 1.722337e-05
hampel   0.9473817 -0.04105608 3.291664e-05
> 
> ## Nice plots -- and check derivatives ----
> 
> head(x. <- seq(-5, 10, length=1501))
[1] -5.00 -4.99 -4.98 -4.97 -4.96 -4.95
> ## [separate lines, for interactive "play": ]
> stopifnot(chkPsiDeriv(p.psiFun(x., "LQQ", par=c(-.5,1.5,.95,NA))))
> stopifnot(chkPsiDeriv(p.psiFun(x., "GGW", par= cG)))
> stopifnot(chkPsiDeriv(p.psiFun(x., "optimal", par=2)))
> stopifnot(chkPsiDeriv(p.psiFun(x., "Hampel",
+                                par = ## Default, but rounded:
+                                round(c(1.5, 3.5, 8) * 0.9016085, 1)),
+                       tol = 1e-3))
> 
> stopifnot(chkPsiDeriv(p.psiFun(x., "biweight", par = 4)))
> stopifnot(chkPsiDeriv(p.psiFun(x., "Welsh", par = 1.5)))
> stopifnot(chkPsiDeriv(p.psiFun(x., "huber", par = 1.5),
+                       tol = c(1e-10, 5e-3)))
> ## "huber"-rho via  Mpsi(*, deriv=-1)  was badly wrong till 2018-06
> 
> ## The same 6, all in one plot:
> op <- par(mfrow=c(3,2), mgp = c(1.5, .6, 0), mar = .1+c(3,3,2,.5))
> p.psiFun2(x., "LQQ", par=c(-.5,1.5,.95,NA))
> p.psiFun2(x., "GGW", par= cG)
> p.psiFun2(x., "optimal", par=1.3)
> p.psiFun2(x., "Hampel", par = round(c(1.5, 3.5, 8) * 0.9016085, 1))
> p.psiFun2(x., "biweight", par = 4)
> p.psiFun2(x., "Welsh", par = 1.5)
> par(op)
> 
> 
> ### (2) Test them as  arguments of  lmrob() or  lmrob.control(): -----
> 
> data(aircraft)
> 
> set.seed(1)
> summary(mp0 <- lmrob(Y ~ ., data = aircraft, psi = 'bisquare', method = 'SMDM'))

Call:
lmrob(formula = Y ~ ., data = aircraft, method = "SMDM", psi = "bisquare")
 \--> method = "SMDM"
Residuals:
    Min      1Q  Median      3Q     Max 
-8.5552 -1.8395 -0.2113  2.8205 46.6311 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.8785690  6.5321336   1.053 0.306256    
X1          -3.2192206  1.0907887  -2.951 0.008543 ** 
X2           1.5876658  0.7442079   2.133 0.046912 *  
X3           0.0018266  0.0004293   4.255 0.000477 ***
X4          -0.0008677  0.0003685  -2.355 0.030083 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Robust residual standard error: 5.335 
Multiple R-squared:  0.7958,	Adjusted R-squared:  0.7504 
Convergence in 22 IRWLS iterations

Robustness weights: 
 observation 22 is an outlier with |weight| = 0 ( < 0.0043); 
 3 weights are ~= 1. The remaining 19 ones are summarized as
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.3958  0.8772  0.9738  0.9139  0.9892  0.9972 
Algorithmic parameters: 
       tuning.chi                bb        tuning.psi        refine.tol 
        1.548e+00         5.000e-01         4.685e+00         1.000e-07 
          rel.tol         scale.tol         solve.tol          zero.tol 
        1.000e-07         1.000e-10         1.000e-07         1.000e-10 
      eps.outlier             eps.x warn.limit.reject warn.limit.meanrw 
        4.348e-03         8.399e-08         5.000e-01         5.000e-01 
     nResample         max.it       best.r.s       k.fast.s          k.max 
           500             50              2              1            200 
   maxit.scale      trace.lev            mts     compute.rd      numpoints 
           200              0           1000              0             10 
fast.s.large.n 
          2000 
                  psi           subsampling                   cov 
           "bisquare"         "nonsingular"             ".vcov.w" 
compute.outlier.stats 
               "SMDM" 
seed : int(0) 
> 
> set.seed(2)
> summary(mp1 <- update(mp0, psi = 'optimal'))

Call:
lmrob(formula = Y ~ ., data = aircraft, method = "SMDM", psi = "optimal")
 \--> method = "SMDM"
Residuals:
    Min      1Q  Median      3Q     Max 
-6.6691 -2.4291  0.2249  3.8876 54.2841 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)   
(Intercept)  9.5007403  5.5576768   1.709  0.10455   
X1          -3.0487969  0.9158751  -3.329  0.00374 **
X2           1.2100330  0.6469186   1.870  0.07777 . 
X3           0.0013810  0.0003910   3.532  0.00238 **
X4          -0.0005549  0.0003269  -1.697  0.10687   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Robust residual standard error: 4.545 
Multiple R-squared:  0.8159,	Adjusted R-squared:  0.775 
Convergence in 1 IRWLS iterations

Robustness weights: 
 2 observations c(16,22) are outliers with |weight| = 0 ( < 0.0043); 
 21 weights are ~= 1.
Algorithmic parameters: 
       tuning.chi                bb        tuning.psi        refine.tol 
        4.047e-01         5.000e-01         1.060e+00         1.000e-07 
          rel.tol         scale.tol         solve.tol          zero.tol 
        1.000e-07         1.000e-10         1.000e-07         1.000e-10 
      eps.outlier             eps.x warn.limit.reject warn.limit.meanrw 
        4.348e-03         8.399e-08         5.000e-01         5.000e-01 
     nResample         max.it       best.r.s       k.fast.s          k.max 
           500             50              2              1            200 
   maxit.scale      trace.lev            mts     compute.rd      numpoints 
           200              0           1000              0             10 
fast.s.large.n 
          2000 
                  psi           subsampling                   cov 
            "optimal"         "nonsingular"             ".vcov.w" 
compute.outlier.stats 
               "SMDM" 
seed : int(0) 
> 
> set.seed(3)
> summary(mp2 <- update(mp0, psi = 'ggw'))

Call:
lmrob(formula = Y ~ ., data = aircraft, method = "SMDM", psi = "ggw")
 \--> method = "SMDM"
Residuals:
    Min      1Q  Median      3Q     Max 
-8.4418 -1.7993 -0.1711  2.8466 47.0906 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.8192703  6.5041383   1.048  0.30831    
X1          -3.1718079  1.0869559  -2.918  0.00918 ** 
X2           1.5705706  0.7510236   2.091  0.05096 .  
X3           0.0017983  0.0004300   4.182  0.00056 ***
X4          -0.0008434  0.0003691  -2.285  0.03466 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Robust residual standard error: 5.367 
Multiple R-squared:  0.7942,	Adjusted R-squared:  0.7484 
Convergence in 20 IRWLS iterations

Robustness weights: 
 observation 22 is an outlier with |weight| <= 0.00044 ( < 0.0043); 
 16 weights are ~= 1. The remaining 6 ones are
     3      4     12     16     17     19 
0.9892 0.9891 0.8770 0.4139 0.9796 0.9839 
Algorithmic parameters: 
      tuning.chi1       tuning.chi2       tuning.chi3       tuning.chi4 
       -5.000e-01         1.500e+00                NA         5.000e-01 
               bb       tuning.psi1       tuning.psi2       tuning.psi3 
        5.000e-01        -5.000e-01         1.500e+00         9.500e-01 
      tuning.psi4        refine.tol           rel.tol         scale.tol 
               NA         1.000e-07         1.000e-07         1.000e-10 
        solve.tol          zero.tol       eps.outlier             eps.x 
        1.000e-07         1.000e-10         4.348e-03         8.399e-08 
warn.limit.reject warn.limit.meanrw 
        5.000e-01         5.000e-01 
     nResample         max.it       best.r.s       k.fast.s          k.max 
           500             50              2              1            200 
   maxit.scale      trace.lev            mts     compute.rd      numpoints 
           200              0           1000              0             10 
fast.s.large.n 
          2000 
                  psi           subsampling                   cov 
                "ggw"         "nonsingular"             ".vcov.w" 
compute.outlier.stats 
               "SMDM" 
seed : int(0) 
> 
> set.seed(4)
> summary(mp3 <- update(mp0, psi = 'welsh'))

Call:
lmrob(formula = Y ~ ., data = aircraft, method = "SMDM", psi = "welsh")
 \--> method = "SMDM"
Residuals:
    Min      1Q  Median      3Q     Max 
-8.7243 -1.9199 -0.2471  2.8060 45.9435 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.6404575  6.5552085   1.013 0.324482    
X1          -3.2329194  1.0954988  -2.951 0.008546 ** 
X2           1.6174887  0.7443222   2.173 0.043367 *  
X3           0.0018656  0.0004279   4.360 0.000378 ***
X4          -0.0008941  0.0003680  -2.430 0.025803 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Robust residual standard error: 5.408 
Multiple R-squared:  0.7958,	Adjusted R-squared:  0.7504 
Convergence in 18 IRWLS iterations

Robustness weights: 
 observation 22 is an outlier with |weight| <= 0.0003 ( < 0.0043); 
 2 weights are ~= 1. The remaining 20 ones are summarized as
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.4284  0.8583  0.9701  0.9112  0.9874  0.9985 
Algorithmic parameters: 
       tuning.chi                bb        tuning.psi        refine.tol 
        5.774e-01         5.000e-01         2.110e+00         1.000e-07 
          rel.tol         scale.tol         solve.tol          zero.tol 
        1.000e-07         1.000e-10         1.000e-07         1.000e-10 
      eps.outlier             eps.x warn.limit.reject warn.limit.meanrw 
        4.348e-03         8.399e-08         5.000e-01         5.000e-01 
     nResample         max.it       best.r.s       k.fast.s          k.max 
           500             50              2              1            200 
   maxit.scale      trace.lev            mts     compute.rd      numpoints 
           200              0           1000              0             10 
fast.s.large.n 
          2000 
                  psi           subsampling                   cov 
              "welsh"         "nonsingular"             ".vcov.w" 
compute.outlier.stats 
               "SMDM" 
seed : int(0) 
> 
> set.seed(5)
> summary(mp4 <- update(mp0, psi = 'ggw', tuning.psi = c(-.5, 1.5, 0.85, NA),
+                       tuning.chi = c(-0.5, 1.5, NA, 0.5)))

Call:
lmrob(formula = Y ~ ., data = aircraft, method = "SMDM", psi = "ggw", tuning.psi = c(-0.5, 
    1.5, 0.85, NA), tuning.chi = c(-0.5, 1.5, NA, 0.5))
 \--> method = "SMDM"
Residuals:
    Min      1Q  Median      3Q     Max 
-7.2207 -2.2226  0.3446  3.5745 52.2885 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)   
(Intercept)  9.6540916  6.6414237   1.454  0.16327   
X1          -3.2353135  1.0947329  -2.955  0.00847 **
X2           1.3343505  0.7636515   1.747  0.09762 . 
X3           0.0015256  0.0004619   3.303  0.00395 **
X4          -0.0006913  0.0003903  -1.771  0.09343 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Robust residual standard error: 5.117 
Multiple R-squared:  0.7832,	Adjusted R-squared:  0.7351 
Convergence in 15 IRWLS iterations

Robustness weights: 
 observation 22 is an outlier with |weight| <= 2.8e-08 ( < 0.0043); 
 15 weights are ~= 1. The remaining 7 ones are
      3       4      12      16      17      19      23 
0.87262 0.79602 0.73029 0.06024 0.96761 0.73117 0.97769 
Algorithmic parameters: 
      tuning.chi1       tuning.chi2       tuning.chi3       tuning.chi4 
       -5.000e-01         1.500e+00                NA         5.000e-01 
               bb       tuning.psi1       tuning.psi2       tuning.psi3 
        5.000e-01        -5.000e-01         1.500e+00         8.500e-01 
      tuning.psi4        refine.tol           rel.tol         scale.tol 
               NA         1.000e-07         1.000e-07         1.000e-10 
        solve.tol          zero.tol       eps.outlier             eps.x 
        1.000e-07         1.000e-10         4.348e-03         8.399e-08 
warn.limit.reject warn.limit.meanrw 
        5.000e-01         5.000e-01 
     nResample         max.it       best.r.s       k.fast.s          k.max 
           500             50              2              1            200 
   maxit.scale      trace.lev            mts     compute.rd      numpoints 
           200              0           1000              0             10 
fast.s.large.n 
          2000 
                  psi           subsampling                   cov 
                "ggw"         "nonsingular"             ".vcov.w" 
compute.outlier.stats 
               "SMDM" 
seed : int(0) 
> 
> set.seed(6)
> summary(mp5 <- update(mp0, psi = 'ggw',
+                       tuning.psi = c(-0.5, 1.0, 0.95, NA),
+                       tuning.chi = c(-0.5, 1.0, NA, 0.5)))

Call:
lmrob(formula = Y ~ ., data = aircraft, method = "SMDM", psi = "ggw", tuning.psi = c(-0.5, 
    1, 0.95, NA), tuning.chi = c(-0.5, 1, NA, 0.5))
 \--> method = "SMDM"
Residuals:
    Min      1Q  Median      3Q     Max 
-8.4182 -1.7447 -0.1322  2.8735 47.0376 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.7557666  6.0919101   1.109 0.282039    
X1          -3.1767976  1.0196958  -3.115 0.005974 ** 
X2           1.5756461  0.7050185   2.235 0.038339 *  
X3           0.0018004  0.0004003   4.497 0.000279 ***
X4          -0.0008432  0.0003446  -2.447 0.024897 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Robust residual standard error: 5.502 
Multiple R-squared:  0.7941,	Adjusted R-squared:  0.7484 
Convergence in 19 IRWLS iterations

Robustness weights: 
 21 weights are ~= 1. The remaining 2 ones are
      16       22 
0.423706 0.005042 
Algorithmic parameters: 
      tuning.chi1       tuning.chi2       tuning.chi3       tuning.chi4 
       -5.000e-01         1.000e+00                NA         5.000e-01 
               bb       tuning.psi1       tuning.psi2       tuning.psi3 
        5.000e-01        -5.000e-01         1.000e+00         9.500e-01 
      tuning.psi4        refine.tol           rel.tol         scale.tol 
               NA         1.000e-07         1.000e-07         1.000e-10 
        solve.tol          zero.tol       eps.outlier             eps.x 
        1.000e-07         1.000e-10         4.348e-03         8.399e-08 
warn.limit.reject warn.limit.meanrw 
        5.000e-01         5.000e-01 
     nResample         max.it       best.r.s       k.fast.s          k.max 
           500             50              2              1            200 
   maxit.scale      trace.lev            mts     compute.rd      numpoints 
           200              0           1000              0             10 
fast.s.large.n 
          2000 
                  psi           subsampling                   cov 
                "ggw"         "nonsingular"             ".vcov.w" 
compute.outlier.stats 
               "SMDM" 
seed : int(0) 
> 
> set.seed(7)
> summary(mp6 <- update(mp0, psi = 'hampel'))

Call:
lmrob(formula = Y ~ ., data = aircraft, method = "SMDM", psi = "hampel")
 \--> method = "SMDM"
Residuals:
   Min     1Q Median     3Q    Max 
-8.706 -1.937 -0.234  2.825 46.037 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.4297956  6.7818616   0.948  0.35564    
X1          -3.1885813  1.1366401  -2.805  0.01170 *  
X2           1.6224243  0.7839018   2.070  0.05315 .  
X3           0.0018590  0.0004445   4.182  0.00056 ***
X4          -0.0008851  0.0003832  -2.310  0.03295 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Robust residual standard error: 5.475 
Multiple R-squared:  0.7946,	Adjusted R-squared:  0.7489 
Convergence in 11 IRWLS iterations

Robustness weights: 
 observation 22 is an outlier with |weight| = 0 ( < 0.0043); 
 20 weights are ~= 1. The remaining 2 ones are
    12     16 
0.8504 0.4975 
Algorithmic parameters: 
      tuning.chi1       tuning.chi2       tuning.chi3                bb 
        3.179e-01         7.417e-01         1.695e+00         5.000e-01 
      tuning.psi1       tuning.psi2       tuning.psi3        refine.tol 
        1.352e+00         3.156e+00         7.213e+00         1.000e-07 
          rel.tol         scale.tol         solve.tol          zero.tol 
        1.000e-07         1.000e-10         1.000e-07         1.000e-10 
      eps.outlier             eps.x warn.limit.reject warn.limit.meanrw 
        4.348e-03         8.399e-08         5.000e-01         5.000e-01 
     nResample         max.it       best.r.s       k.fast.s          k.max 
           500             50              2              1            200 
   maxit.scale      trace.lev            mts     compute.rd      numpoints 
           200              0           1000              0             10 
fast.s.large.n 
          2000 
                  psi           subsampling                   cov 
             "hampel"         "nonsingular"             ".vcov.w" 
compute.outlier.stats 
               "SMDM" 
seed : int(0) 
> 
> set.seed(8)
> ctr7 <- lmrob.control(psi = 'ggw',
+                       tuning.psi = c(-0.3, 1.4, 0.95, NA),
+                       tuning.chi = c(-0.3, 1.4, NA, 0.5))
> ctr7$tuning.psi ## -> "constants"
[1] -0.30  1.40  0.95    NA
attr(,"constants")
[1] 0.0000000 2.0011562 1.4000000 0.4125717 5.6874488
> ctr7$tuning.chi
[1] -0.3  1.4   NA  0.5
attr(,"constants")
[1] 0.00000000 0.24044569 1.40000000 0.09081713 0.27558437
> summary(mp7 <-lmrob(Y ~ ., data = aircraft, control = ctr7)) # *not* converging in k.max=200

Call:
lmrob(formula = Y ~ ., data = aircraft, control = ctr7)
 \--> method = "S"
Residuals:
    Min      1Q  Median      3Q     Max 
-7.6919 -1.9269  0.1767  3.7081 48.5801 

Algorithm did not converge

Coefficients of the *initial* S-estimator:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 13.155499         NA      NA       NA
X1          -4.349383         NA      NA       NA
X2           1.647243         NA      NA       NA
X3           0.001817         NA      NA       NA
X4          -0.001035         NA      NA       NA

Robustness weights: 
 2 observations c(16,22) are outliers with |weight| <= 0.0003 ( < 0.0043); 
 4 weights are ~= 1. The remaining 17 ones are summarized as
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.03668 0.20120 0.58420 0.52290 0.71930 0.99110 
Algorithmic parameters: 
      tuning.chi1       tuning.chi2       tuning.chi3       tuning.chi4 
       -3.000e-01         1.400e+00                NA         5.000e-01 
               bb       tuning.psi1       tuning.psi2       tuning.psi3 
        5.000e-01        -3.000e-01         1.400e+00         9.500e-01 
      tuning.psi4        refine.tol           rel.tol         scale.tol 
               NA         1.000e-07         1.000e-07         1.000e-10 
        solve.tol          zero.tol       eps.outlier             eps.x 
        1.000e-07         1.000e-10         4.348e-03         8.399e-08 
warn.limit.reject warn.limit.meanrw 
        5.000e-01         5.000e-01 
     nResample         max.it       best.r.s       k.fast.s          k.max 
           500             50              2              1            200 
   maxit.scale      trace.lev            mts     compute.rd fast.s.large.n 
           200              0           1000              0           2000 
                  psi           subsampling                   cov 
                "ggw"         "nonsingular"         ".vcov.avar1" 
compute.outlier.stats 
                 "SM" 
seed : int(0) 
Warning messages:
1: In lmrob.S(x, y, control = control) :
  S refinements did not converge (to refine.tol=1e-07) in 200 (= k.max) steps
2: In lmrob.fit(x, y, control, init = init) :
  initial estim. 'init' not converged -- will be return()ed basically unchanged
> 
> set.seed(9)
> summary(mp8 <- update(mp0, psi = 'lqq'))

Call:
lmrob(formula = Y ~ ., data = aircraft, method = "SMDM", psi = "lqq")
 \--> method = "SMDM"
Residuals:
   Min     1Q Median     3Q    Max 
-8.280 -1.717 -0.138  2.857 47.743 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  7.0858566  6.3506262   1.116 0.279194    
X1          -3.1657682  1.0600204  -2.987 0.007914 ** 
X2           1.5402736  0.7336570   2.099 0.050145 .  
X3           0.0017612  0.0004222   4.171 0.000574 ***
X4          -0.0008188  0.0003616  -2.265 0.036118 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Robust residual standard error: 5.319 
Multiple R-squared:  0.7944,	Adjusted R-squared:  0.7487 
Convergence in 19 IRWLS iterations

Robustness weights: 
 observation 22 is an outlier with |weight| = 0 ( < 0.0043); 
 16 weights are ~= 1. The remaining 6 ones are
     3      4     12     16     17     19 
0.9861 0.9842 0.8921 0.3720 0.9820 0.9782 
Algorithmic parameters: 
      tuning.chi1       tuning.chi2       tuning.chi3       tuning.chi4 
       -5.000e-01         1.500e+00                NA         5.000e-01 
               bb       tuning.psi1       tuning.psi2       tuning.psi3 
        5.000e-01        -5.000e-01         1.500e+00         9.500e-01 
      tuning.psi4        refine.tol           rel.tol         scale.tol 
               NA         1.000e-07         1.000e-07         1.000e-10 
        solve.tol          zero.tol       eps.outlier             eps.x 
        1.000e-07         1.000e-10         4.348e-03         8.399e-08 
warn.limit.reject warn.limit.meanrw 
        5.000e-01         5.000e-01 
     nResample         max.it       best.r.s       k.fast.s          k.max 
           500             50              2              1            200 
   maxit.scale      trace.lev            mts     compute.rd      numpoints 
           200              0           1000              0             10 
fast.s.large.n 
          2000 
                  psi           subsampling                   cov 
                "lqq"         "nonsingular"             ".vcov.w" 
compute.outlier.stats 
               "SMDM" 
seed : int(0) 
> 
> set.seed(10) ##  c(.) drops attributes :
> ctr9 <- lmrob.control(psi = 'lqq', tuning.psi = c(ctr7$tuning.psi), tuning.chi = c(ctr7$tuning.chi))
> ctr9$tuning.psi
[1] -0.30  1.40  0.95    NA
attr(,"constants")
[1] 1.3007171 0.9290836 1.3000000
> ctr9$tuning.chi
[1] -0.3  1.4   NA  0.5
attr(,"constants")
[1] 0.2763568 0.1973977 1.3000000
> ## Confirm these constants above (against the ones we got earlier)
> ## by recomputing them using higher accuracy :
> (tpsi. <- do.call(.psi.lqq.findc, c(ctr9$tuning.psi, list(rel.tol=1e-11, tol=1e-8))))
[1] 1.3007495 0.9291068 1.3000000
> (tchi. <- do.call(.psi.lqq.findc, c(ctr9$tuning.chi, list(rel.tol=1e-11, tol=1e-8))))
[1] 0.2763425 0.1973875 1.3000000
> (tol4 <- .Machine$double.eps^.25)
[1] 0.0001220703
> 
> Rver <- getRversion()
> integr.bug <- "2.12.0" <= Rver && Rver <= "3.0.1"
> integr.bug
[1] FALSE
> if(integr.bug) tol4 <- 8*tol4
> 
> assert.EQ(attr(ctr9$tuning.psi, "constants"), tpsi., tol=tol4, giveRE=TRUE)
Mean relative difference: 2.495013e-05
> assert.EQ(attr(ctr9$tuning.chi, "constants"), tchi., tol=tol4, giveRE=TRUE)
Mean relative difference: 5.155651e-05
> 
> summary(mp9 <- lmrob(Y ~ ., data = aircraft, control = ctr9))

Call:
lmrob(formula = Y ~ ., data = aircraft, control = ctr9)
 \--> method = "MM"
Residuals:
     Min       1Q   Median       3Q      Max 
-10.4061  -2.6517  -0.4156   3.7945  38.6444 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)   
(Intercept)  3.5003005 12.9625202   0.270  0.79021   
X1          -3.2953770  0.9467913  -3.481  0.00267 **
X2           1.8957842  0.9928099   1.910  0.07227 . 
X3           0.0022793  0.0014340   1.589  0.12936   
X4          -0.0011563  0.0008966  -1.290  0.21347   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Robust residual standard error: 6.127 
Multiple R-squared:  0.7973,	Adjusted R-squared:  0.7523 
Convergence in 33 IRWLS iterations

Robustness weights: 
 17 weights are ~= 1. The remaining 6 ones are
      3       4      12      16      17      22 
0.97698 0.99840 0.82584 0.78662 0.91318 0.06838 
Algorithmic parameters: 
      tuning.chi1       tuning.chi2       tuning.chi3       tuning.chi4 
       -3.000e-01         1.400e+00                NA         5.000e-01 
               bb       tuning.psi1       tuning.psi2       tuning.psi3 
        5.000e-01        -3.000e-01         1.400e+00         9.500e-01 
      tuning.psi4        refine.tol           rel.tol         scale.tol 
               NA         1.000e-07         1.000e-07         1.000e-10 
        solve.tol          zero.tol       eps.outlier             eps.x 
        1.000e-07         1.000e-10         4.348e-03         8.399e-08 
warn.limit.reject warn.limit.meanrw 
        5.000e-01         5.000e-01 
     nResample         max.it       best.r.s       k.fast.s          k.max 
           500             50              2              1            200 
   maxit.scale      trace.lev            mts     compute.rd fast.s.large.n 
           200              0           1000              0           2000 
                  psi           subsampling                   cov 
                "lqq"         "nonsingular"         ".vcov.avar1" 
compute.outlier.stats 
                 "SM" 
seed : int(0) 
> 
> 
> cat('Time elapsed: ', proc.time(),'\n') # for ``statistical reasons''
Time elapsed:  3.35 0.28 3.62 NA NA 
> 
> proc.time()
   user  system elapsed 
   3.35    0.28    3.62