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) > > L1.obj <- function(m, X) sum(sqrt(rowSums((X - rep(m, each = nrow(X)))^2))) > > ## Really 'simple' 3-D example > X <- rbind(cbind(1:12, round(abs(-5:6)^1.5) %% 7, rep(1:4,3)), + c(3,0, 20), c(6,1, 200))## 'outliers' > > c0 <- system.time(for(i in 1:10) m0 <- L1median(X, method="nlminb",tol =3e-16)) > c1 <- system.time(for(i in 1:10) m1 <- L1median(X, method="nlm", tol = 3e-16)) > c2 <- system.time(for(i in 1:10) m2 <- L1median(X, method="HoCr", tol = 3e-16)) > c3 <- system.time(for(i in 1:10) m3 <- L1median(X, method="Vardi",tol = 3e-16)) > (cM <- rbind(nlminb = c0, nlm = c1, HoCr = c2, Vardi = c3)) user.self sys.self elapsed user.child sys.child nlminb 0.03 0 0.03 NA NA nlm 0.01 0 0.01 NA NA HoCr 0.02 0 0.01 NA NA Vardi 0.02 0 0.02 NA NA > mNms <- rownames(cM) # method names > > pX <- matrix(NA, ncol(X), length(mNms), + dimnames = list(NULL, mNms)) > ## FIXME, once we have common return value > pX[,"nlminb"] <- m0$par > pX[,"nlm" ] <- m1$estimate > pX[,"HoCr" ] <- m2 > pX[,"Vardi" ] <- m3 > t(pX) [,1] [,2] [,3] nlminb 6.563449 2.017194 3.164644 nlm 6.563445 2.017189 3.164644 HoCr 6.563445 2.017189 3.164644 Vardi 6.563445 2.017189 3.164644 > > apply(pX, 2, function(m) L1.obj(m,X)) - 259.7393299943 nlminb nlm HoCr Vardi 9.595169e-11 4.359890e-11 4.359890e-11 4.359890e-11 > ## hmm, "nlminb" is slightly "worse" than the others .. > > for(j in 1:ncol(pX)) + stopifnot(all.equal(pX[,j], pX[,1 + j %% ncol(pX)], tol = 1e-6)) > ## hmm: 1e-6 currently needed 'cause of "nlminb" > > ##--- even much simpler: Examples where coord.wise.median == data point > (x3 <- cbind(c(0,1,8),0:2)) [,1] [,2] [1,] 0 0 [2,] 1 1 [3,] 8 2 > x5 <- rbind(x3, + c(1,0), + c(2,1)) > (x5 <- x5[order(x5[,1],x5[,2]), ]) [,1] [,2] [1,] 0 0 [2,] 1 0 [3,] 1 1 [4,] 2 1 [5,] 8 2 > m3.0 <- L1median(x3, method="nlminb", trace=2, tol=3e-16) pscale= 1,1 0: 8.4852814: 1.00000 1.00000 > m3.1 <- L1median(x3, method="nlm", trace=2, tol=3e-16) pscale= 1,1 iteration = 0 Parameter: [1] 1 1 Function Value [1] 8.485281 Gradient: [1] NaN NaN Relative gradient close to zero. Current iterate is probably solution. > m3.2 <- L1median(x3, method="HoCrJo", trace=2) k= 1 obj= 8.48528137526 m=(1,1) nstep=28 halvings > m3.3 <- L1median(x3, method="Vardi", trace=2) pscale= 1,1 y == x_2 needed 1 iterations (0 of which had y==x_k) > > pX <- matrix(NA, ncol(x3), length(mNms), dimnames = list(NULL, mNms)) > pX[,"nlminb"] <- m3.0$par > pX[,"nlm" ] <- m3.1$estimate > pX[,"HoCr" ] <- m3.2 > pX[,"Vardi" ] <- m3.3 > t(pX) [,1] [,2] nlminb 1 1 nlm 1 1 HoCr 1 1 Vardi 1 1 > for(j in 1:ncol(pX)) + stopifnot(all.equal(pX[,j], pX[,1 + j %% ncol(pX)], tol = 1e-6)) > > > m5.0 <- L1median(x5, method="nlminb", trace=2, tol=3e-16) pscale= 1,1 0: 10.485281: 1.00000 1.00000 2: 10.252811: 1.26833 0.593349 4: 10.245887: 1.27176 0.664737 6: 10.245883: 1.27004 0.664946 > m5.1 <- L1median(x5, method="nlm", trace=2, tol=3e-16) # fails !! pscale= 1,1 iteration = 0 Parameter: [1] 1 1 Function Value [1] 10.48528 Gradient: [1] NaN NaN Relative gradient close to zero. Current iterate is probably solution. > m5.2 <- L1median(x5, method="HoCrJo", trace=2) k= 1 obj= 10.3739361299 m=(1.45,0.4504) k= 2 obj= 10.260048966 m=(1.329,0.5898) k= 3 obj= 10.2484019202 m=(1.293,0.6316) k= 4 obj= 10.246407523 m=(1.28,0.6491) k= 5 obj= 10.2459999303 m=(1.274,0.6572) k= 6 obj= 10.2459097722 m=(1.272,0.6611) k= 7 obj= 10.245889006 m=(1.271,0.663) k= 8 obj= 10.2458841074 m=(1.27,0.664) k= 9 obj= 10.2458829325 m=(1.27,0.6645) k= 10 obj= 10.2458826467 m=(1.27,0.6647) k= 11 obj= 10.2458825764 m=(1.27,0.6648) k= 12 obj= 10.2458825588 m=(1.27,0.6649) k= 13 obj= 10.2458825544 m=(1.27,0.6649) k= 14 obj= 10.2458825533 m=(1.27,0.6649) k= 15 obj= 10.245882553 m=(1.27,0.6649) k= 16 obj= 10.2458825529 m=(1.27,0.6649) k= 17 obj= 10.2458825529 m=(1.27,0.6649) k= 18 obj= 10.2458825529 m=(1.27,0.6649) k= 19 obj= 10.2458825529 m=(1.27,0.6649) k= 20 obj= 10.2458825529 m=(1.27,0.6649) k= 21 obj= 10.2458825529 m=(1.27,0.6649) k= 22 obj= 10.2458825529 m=(1.27,0.6649) k= 23 obj= 10.2458825529 m=(1.27,0.6649) k= 24 obj= 10.2458825529 m=(1.27,0.6649) nstep= 3 halvings > m5.3 <- L1median(x5, method="Vardi", trace=2) pscale= 1,1 y == x_3 k= 1 rel.chg= 0.2594973, m=(1.2279, 0.7219) k= 2 rel.chg= 0.02412694, m=(1.2484, 0.69554) k= 3 rel.chg= 0.01337325, m=(1.2595, 0.68074) k= 4 rel.chg= 0.006903148, m=(1.2651, 0.67294) k= 5 rel.chg= 0.003431709, m=(1.2678, 0.66896) k= 6 rel.chg= 0.001674532, m=( 1.269, 0.66696) k= 7 rel.chg=0.0008101124, m=(1.2696, 0.66595) k= 8 rel.chg=0.0003905495, m=(1.2698, 0.66545) k= 9 rel.chg=0.0001881044, m=( 1.27, 0.6652) k= 10 rel.chg= 9.06322e-05, m=( 1.27, 0.66507) k= 11 rel.chg=4.371475e-05, m=( 1.27, 0.66501) k= 12 rel.chg= 2.11158e-05, m=( 1.27, 0.66498) k= 13 rel.chg=1.021718e-05, m=( 1.27, 0.66496) k= 14 rel.chg=4.953111e-06, m=( 1.27, 0.66495) k= 15 rel.chg=2.406099e-06, m=( 1.27, 0.66495) k= 16 rel.chg= 1.17137e-06, m=( 1.27, 0.66495) k= 17 rel.chg=5.715738e-07, m=( 1.27, 0.66494) k= 18 rel.chg=3.043732e-07, m=( 1.27, 0.66494) k= 19 rel.chg=1.629409e-07, m=( 1.27, 0.66494) k= 20 rel.chg=8.706675e-08, m=( 1.27, 0.66494) k= 21 rel.chg=4.644881e-08, m=( 1.27, 0.66494) k= 22 rel.chg=2.474479e-08, m=( 1.27, 0.66494) k= 23 rel.chg=1.316601e-08, m=( 1.27, 0.66494) needed 24 iterations (0 of which had y==x_k) > > pX <- matrix(NA, ncol(x5), length(mNms), dimnames = list(NULL, mNms)) > pX[,"nlminb"] <- m5.0$par > pX[,"nlm" ] <- m5.1$estimate > pX[,"HoCr" ] <- m5.2 > pX[,"Vardi" ] <- m5.3 > t(pX) [,1] [,2] nlminb 1.270045 0.6649437 nlm 1.000000 1.0000000 HoCr 1.270045 0.6649437 Vardi 1.270045 0.6649437 > > ## FIXME; nlm() currently fails here: > pX <- pX[, - which("nlm" == mNms)] > > for(j in 1:ncol(pX)) + stopifnot(all.equal(pX[,j], pX[,1 + j %% ncol(pX)], tol = 1e-6)) > > > > stopifnot(require(MASS)) Loading required package: MASS > ##-> lazy loading data sets > > (sl.HC <- L1median (stackloss, trace = TRUE, method = "HoCrJo")) needed 26 iterations with a total of 4 stepsize halvings Air.Flow Water.Temp Acid.Conc. stack.loss 59.03170 20.68484 86.66082 15.51665 > > system.time(rr0 <- L1median(stackloss, method="HoCrJo", tol = 1e-14)) user system elapsed 0.02 0.00 0.02 > system.time(rr1 <- L1median(stackloss, method="Nelder-M", tol = 1e-14)) user system elapsed 0 0 0 > > system.time(rr2 <- L1median(stackloss, method="BFGS", tol = 1e-14)) user system elapsed 0 0 0 > ## MM: Hmm, this ("CG") now takes MUCH longer (factor 5 - 10 !): > system.time(rr3 <- L1median(stackloss, method="CG", tol = 1e-14)) user system elapsed 0.02 0.00 0.02 > ## takes even longer: > system.time(rr3.2 <- L1median(stackloss, method="CG", tol = 1e-14, type = 2)) user system elapsed 0 0 0 > ## (fastest): > system.time(rr3.3 <- L1median(stackloss, method="CG", tol = 1e-14, type = 3)) user system elapsed 0 0 0 > > ## nlm with gradient: > system.time(rr4 <- L1median(stackloss, method="nlm", tol = 1e-16)) user system elapsed 0 0 0 > ##--> fastest! {faster than rr0 by almost factor of two!} > system.time(rr4 <- L1median(stackloss, method="nlm", tol = 1e-16, trace = 2)) pscale= 4.8,2.2,3.2,4.4 iteration = 0 Step: [1] 0 0 0 0 Parameter: [1] 58 20 87 15 Function Value [1] 255.2048 Gradient: [1] -1.9352118 -0.9406718 0.9142917 -0.3972367 iteration = 1 Step: [1] 1.0624342 0.1084862 -0.2230878 0.1832507 Parameter: [1] 59.06243 20.10849 86.77691 15.18325 Function Value [1] 254.0452 Gradient: [1] 0.3262866 -1.0379401 0.2693502 -0.5045123 iteration = 2 Step: [1] -0.04427560 0.12467870 -0.08745203 0.23919501 Parameter: [1] 59.01816 20.23316 86.68946 15.42245 Function Value [1] 253.8375 Gradient: [1] 0.10437905 -0.87309103 0.08669800 -0.02535936 iteration = 3 Step: [1] -0.01906074 0.23737083 -0.09451050 0.13879161 Parameter: [1] 58.99910 20.47054 86.59495 15.56124 Function Value [1] 253.6938 Gradient: [1] -0.0362155 -0.4390493 -0.1179374 0.1928198 iteration = 4 Step: [1] 0.005322241 0.155754249 -0.007711965 -0.012900326 Parameter: [1] 59.00442 20.62629 86.58724 15.54834 Function Value [1] 253.6491 Gradient: [1] -0.04719098 -0.11993896 -0.14283783 0.10873757 iteration = 5 Step: [1] 0.01274755 0.05423190 0.03565277 -0.03514507 Parameter: [1] 59.01717 20.68052 86.62289 15.51319 Function Value [1] 253.6395 Gradient: [1] -0.019724233 -0.002549225 -0.073443635 0.005786863 iteration = 6 Step: [1] 0.009132417 0.005710474 0.025190856 -0.005267043 Parameter: [1] 59.02630 20.68623 86.64808 15.50792 Function Value [1] 253.6382 Gradient: [1] -0.004229506 0.007627714 -0.023994298 -0.014666831 iteration = 7 Step: [1] 0.004674191 -0.000932607 0.011407086 0.005121485 Parameter: [1] 59.03097 20.68530 86.65949 15.51305 Function Value [1] 253.638 Gradient: [1] 0.0005902791 0.0024114084 -0.0022162298 -0.0072083009 iteration = 8 Step: [1] 0.0007631427 -0.0004833079 0.0014954058 0.0030731835 Parameter: [1] 59.03174 20.68482 86.66098 15.51612 Function Value [1] 253.638 Gradient: [1] 0.0003305839 0.0001253512 0.0003865554 -0.0011391887 iteration = 9 Step: [1] -2.551964e-05 5.638316e-06 -1.325580e-04 5.078933e-04 Parameter: [1] 59.03171 20.68482 86.66085 15.51663 Function Value [1] 253.638 Gradient: [1] 3.561266e-05 -2.959142e-05 6.990797e-05 -4.767458e-05 iteration = 10 Step: [1] -1.324359e-05 1.428300e-05 -3.290767e-05 2.056096e-05 Parameter: [1] 59.03170 20.68484 86.66082 15.51665 Function Value [1] 253.638 Gradient: [1] 1.469642e-06 -3.596113e-06 2.698806e-06 8.180226e-07 iteration = 11 Step: [1] -7.930822e-07 1.511344e-06 -1.399551e-06 -3.631488e-07 Parameter: [1] 59.03170 20.68484 86.66082 15.51665 Function Value [1] 253.638 Gradient: [1] 3.180428e-08 -1.849170e-07 -4.009023e-08 8.478721e-08 iteration = 12 Step: [1] -1.359474e-08 8.257469e-08 1.572360e-08 -2.893809e-08 Parameter: [1] 59.03170 20.68484 86.66082 15.51665 Function Value [1] 253.638 Gradient: [1] 8.071191e-10 -6.169280e-09 -7.256855e-09 5.295721e-10 iteration = 13 Step: [1] 3.910685e-10 3.274071e-09 3.616336e-09 4.987015e-10 Parameter: [1] 59.03170 20.68484 86.66082 15.51665 Function Value [1] 253.638 Gradient: [1] 2.407144e-11 -6.649258e-12 -2.216356e-10 -1.507106e-10 iteration = 14 Step: [1] 2.433609e-11 2.498268e-11 1.154632e-10 8.559908e-11 Parameter: [1] 59.03170 20.68484 86.66082 15.51665 Function Value [1] 253.638 Gradient: [1] 1.194461e-12 2.523044e-12 -1.822556e-12 -4.373057e-12 iteration = 15 Step: [1] -2.842171e-14 -8.633094e-13 9.663381e-13 1.961098e-12 Parameter: [1] 59.03170 20.68484 86.66082 15.51665 Function Value [1] 253.638 Gradient: [1] 5.928591e-14 3.386874e-14 1.262879e-15 -3.490264e-14 iteration = 16 Step: [1] -2.842171e-14 -1.776357e-14 0.000000e+00 7.105427e-15 Parameter: [1] 59.03170 20.68484 86.66082 15.51665 Function Value [1] 253.638 Gradient: [1] -3.608225e-16 6.383782e-16 7.514822e-15 1.311451e-15 iteration = 17 Step: [1] 0.000000e+00 0.000000e+00 0.000000e+00 -1.776357e-15 Parameter: [1] 59.03170 20.68484 86.66082 15.51665 Function Value [1] 253.638 Gradient: [1] 6.106227e-16 1.346145e-15 7.487067e-15 -2.234324e-15 iteration = 18 Parameter: [1] 59.03170 20.68484 86.66082 15.51665 Function Value [1] 253.638 Gradient: [1] 6.106227e-16 1.346145e-15 7.487067e-15 -2.234324e-15 Successive iterates within tolerance. Current iterate is probably solution. user system elapsed 0 0 0 > > mm <- rbind(HoCrJo= rr0, "NelderM"=rr1$par, BFGS=rr2$par, + CG = rr3$par, CG.2 = rr3.2$par, CG.3= rr3.3$par, + "nlm.grad" = rr4$estimate) > L1.obj <- function(m, Xmat) sum(sqrt(rowSums((Xmat - rep(m, each = nrow(Xmat)))^2))) > mm <- cbind(mm, Obj = apply(mm, 1, L1.obj, Xmat = stackloss)) > print(mm[,"Obj"] - min(mm[,"Obj"]), digits = 4) HoCrJo NelderM BFGS CG CG.2 CG.3 nlm.grad 0.000e+00 9.202e-08 1.137e-13 1.705e-13 1.137e-13 8.527e-14 0.000e+00 > ## HoCrJo is best (since it iterates too long!) *together* with nlm; > ## but only significantly wrt "NelderM" > op <- options(digits=12)# more than usual to see (irrelevant) differences: > mm Air.Flow Water.Temp Acid.Conc. stack.loss Obj HoCrJo 59.0316977522 20.6848379013 86.6608169479 15.5166475602 253.637976379 NelderM 59.0317619939 20.6848293601 86.6605411339 15.5165691054 253.637976471 BFGS 59.0316980520 20.6848379843 86.6608169146 15.5166479785 253.637976379 CG 59.0316982171 20.6848380718 86.6608169461 15.5166479708 253.637976379 CG.2 59.0316981085 20.6848380609 86.6608169320 15.5166479564 253.637976379 CG.3 59.0316981155 20.6848379508 86.6608169948 15.5166478824 253.637976379 nlm.grad 59.0316978375 20.6848379675 86.6608169883 15.5166476489 253.637976379 > options(op) > > > proc.time() user system elapsed 0.54 0.09 0.62