R Under development (unstable) (2024-06-18 r86781 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. > stopifnot(require("robustbase")) Loading required package: robustbase > source(system.file("xtraR", "platform-sessionInfo.R", # moreSessionInfo() etc + package = "robustbase", mustWork=TRUE)) > ## testing functions: > source(system.file("xtraR/test-tools.R", package = "robustbase")) > ## showProc.time(), showSys.time() ... > S.time <- showSys.time # "back compatible" > > mS <- moreSessionInfo(print.=TRUE) List of 5 $ sizeof.long : int 4 $ sizeof.longlong : int 8 $ sizeof.longdouble: int 16 $ sizeof.pointer : int 8 $ sizeof.time_t : int 8 64 bit platform type 'windows' ==> onWindows: TRUE arch: x86-64 osVersion (0): Windows Server 2022 x64 (build 20348) osVersion: 'Windows Server 2022 x64 (build 20348)' + BLAS "is" Lapack: TRUE | BLAS=OpenBLAS: FALSE | Lapack=OpenBLAS: FALSE strictR: FALSE > > ## as long as we don't export these (nor provide an nlrob(., method=.) interface: > nlrob.MM <- robustbase:::nlrob.MM > nlrob.tau <- robustbase:::nlrob.tau > nlrob.CM <- robustbase:::nlrob.CM > nlrob.mtl <- robustbase:::nlrob.mtl > > (doExtras <- robustbase:::doExtras()) [1] FALSE > if(doExtras) { + NP <- 30 ; tol <- 1e-11 + } else { ## "fast" + NP <- 15 ; tol <- 1e-7 + } > > start.from.true <- !doExtras # (but not necessarily ..) > if(start.from.true) { # population size = NP (random) + 1 (true parameters) + init_p <- c(1, 0.2) + init_p_sigma <- c(1, 0.2, 1) + } else { + init_p <- init_p_sigma <- NULL + } > > if(!dev.interactive(orNone=TRUE)) pdf("nlregrob-tst.pdf") > > RNGversion("3.5.0") # -- TODO once R >> 3.5.0 : update results !! Warning message: In RNGkind("Mersenne-Twister", "Inversion", "Rounding") : non-uniform 'Rounding' sampler used > > ## Stromberg, Arnold J. (1993). > ## Computation of high breakdown nonlinear regression parameters. > ## J. Amer. Statist. Assoc. 88(421), 237-244. > > ## exponential regression > Expo <- function(x, a, b) exp(a + b*x) > set.seed(2345) # for reproducibility > ## data without outliers: > d.exp30 <- data.frame(x = sort( runif(30, 0, 10) ), err = rnorm(30)) > d.exp30 <- transform(d.exp30, y = Expo(x, 1, 0.2) + err) > ## classical (starting at truth .. hmm) > op <- options(digits=12) > Cfit <- nls(y ~ Expo(x, a, b), data = d.exp30, start = c(a = 1, b = 0.2), + control = nls.control(tol = 8e-8, printEval = TRUE), trace=TRUE) 19.5256274928 (5.14e-02): par = (1 0.2) It. 1, fac= 1, eval (no.,total): ( 1, 1): new dev = 19.4745 19.4745496771 (7.25e-04): par = (0.996314343208 0.201072709229) It. 2, fac= 1, eval (no.,total): ( 1, 2): new dev = 19.4745 19.4745396962 (1.74e-05): par = (0.996513207557 0.201043400501) It. 3, fac= 1, eval (no.,total): ( 1, 3): new dev = 19.4745 19.4745396904 (4.31e-07): par = (0.996508092459 0.201044115571) It. 4, fac= 1, eval (no.,total): ( 1, 4): new dev = 19.4745 19.4745396904 (2.03e-08): par = (0.996508219071 0.201044097876) > showProc.time()# ---- OS X needing 6e-8 Time (user system elapsed): 0.09 0.01 0.11 > options(op) > > ## robust > Rfit.MM.S.bisquare <- + nlrob.MM(y ~ Expo(x, a, b), data = d.exp30, + lower = c(a = -10, b = -2), upper = c(10, 2), + NP = NP, tol = tol, add_to_init_pop = init_p ) > if(doExtras) { + Rfit.MM.S.lqq <- update(Rfit.MM.S.bisquare, psi = "lqq") + Rfit.MM.S.optimal <- update(Rfit.MM.S.bisquare, psi = "optimal") + Rfit.MM.S.hampel <- update(Rfit.MM.S.bisquare, psi = "hampel") + } > showProc.time() Time (user system elapsed): 0.41 0 0.41 > Rfit.MM.lts.bisquare <- update(Rfit.MM.S.bisquare, init = "lts") > Rfit.MM.lts.lqq <- update(Rfit.MM.S.bisquare, init = "lts", psi = "lqq") > Rfit.MM.lts.optimal <- update(Rfit.MM.S.bisquare, init = "lts", psi = "optimal") > Rfit.MM.lts.hampel <- update(Rfit.MM.S.bisquare, init = "lts", psi = "hampel") > showProc.time() Time (user system elapsed): 0.58 0.02 0.59 > > S.time(Rfit.tau.bisquare <- + nlrob.tau( y ~ Expo(x, a, b), data = d.exp30, + lower = c(a = -10, b = -2), upper = c(10, 2), + NP = NP, add_to_init_pop = init_p )) Time user system elapsed Time 0.43 0.00 0.42 > S.time(Rfit.tau.optimal <- update(Rfit.tau.bisquare, psi = "optimal")) Time user system elapsed Time 0.35 0.00 0.34 > > S.time(Rfit.CM <- nlrob.CM( y ~ Expo(x, a, b), data = d.exp30, + lower = c(a = -10, b = -2, sigma = 0), + upper = c( 10, 2, 10), + NP = NP, add_to_init_pop = init_p_sigma )) Time user system elapsed Time 0.32 0.00 0.31 > S.time(Rfit.mtl <- nlrob.mtl(y ~ Expo(x, a, b), data = d.exp30, + lower = c(a = -10, b = -2, sigma = 0), + upper = c( 10, 2, 3), + NP = ceiling(NP*1.33), # <- higher prob. to get close + tol = tol, + trace=TRUE, details=TRUE, + add_to_init_pop = init_p_sigma )) 1 : < 0.1744573 > ( 74.66194 ) 1 0.2 1 2 : < 0.121763 > ( 74.66194 ) 1 0.2 1 3 : < 0.09241423 > ( 74.66194 ) 1 0.2 1 4 : < 0.08943249 > ( 74.66194 ) 1 0.2 1 5 : < 0.07877094 > ( 74.66194 ) 1 0.2 1 6 : < 0.0780325 > ( 74.66194 ) 1 0.2 1 7 : < 0.07270791 > ( 74.66194 ) 1 0.2 1 8 : < 0.07212522 > ( 74.66194 ) 1 0.2 1 9 : < 0.06919738 > ( 64.39819 ) 1.133355 0.09967391 0.9831099 10 : < 0.06455391 > ( 63.93903 ) 1.189944 0.0741118 0.9811784 11 : < 0.06455391 > ( 63.93903 ) 1.189944 0.0741118 0.9811784 12 : < 0.06455391 > ( 63.93903 ) 1.189944 0.0741118 0.9811784 13 : < 0.06455391 > ( 63.93903 ) 1.189944 0.0741118 0.9811784 14 : < 0.06104868 > ( 63.93903 ) 1.189944 0.0741118 0.9811784 15 : < 0.06104868 > ( 63.93903 ) 1.189944 0.0741118 0.9811784 16 : < 0.06104868 > ( 63.93903 ) 1.189944 0.0741118 0.9811784 17 : < 0.04815709 > ( 63.93903 ) 1.189944 0.0741118 0.9811784 18 : < 0.03052925 > ( 63.93903 ) 1.189944 0.0741118 0.9811784 19 : < 0.03129668 > ( 63.49528 ) 1 0.245802 1 20 : < 0.03129668 > ( 63.49528 ) 1 0.245802 1 21 : < 0.03129668 > ( 63.49528 ) 1 0.245802 1 22 : < 0.06921241 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 23 : < 0.06921241 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 24 : < 0.06580972 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 25 : < 0.0656521 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 26 : < 0.06081836 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 27 : < 0.06081836 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 28 : < 0.06081836 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 29 : < 0.06081836 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 30 : < 0.06081836 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 31 : < 0.06081836 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 32 : < 0.06081836 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 33 : < 0.06000373 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 34 : < 0.06000373 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 35 : < 0.06000373 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 36 : < 0.06000373 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 37 : < 0.06000373 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 38 : < 0.06000373 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 39 : < 0.06000373 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 40 : < 0.05956777 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 41 : < 0.04766818 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 42 : < 0.04766818 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 43 : < 0.04766818 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 44 : < 0.04766818 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 45 : < 0.04766818 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 46 : < 0.04766818 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 47 : < 0.04161718 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 48 : < 0.04161718 > ( 41.57127 ) 0.5981937 0.3009903 0.6217691 49 : < 0.05793346 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 50 : < 0.05793346 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 51 : < 0.05793346 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 52 : < 0.05793346 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 53 : < 0.05786467 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 54 : < 0.05786467 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 55 : < 0.05786467 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 56 : < 0.05786467 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 57 : < 0.05786467 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 58 : < 0.05786467 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 59 : < 0.05786467 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 60 : < 0.0531622 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 61 : < 0.0531622 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 62 : < 0.0531622 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 63 : < 0.0531622 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 64 : < 0.0531622 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 65 : < 0.0531622 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 66 : < 0.0531622 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 67 : < 0.05315922 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 68 : < 0.05133014 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 69 : < 0.05133014 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 70 : < 0.05133014 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 71 : < 0.05133014 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 72 : < 0.0503882 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 73 : < 0.05006865 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 74 : < 0.05006865 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 75 : < 0.05006865 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 76 : < 0.04752007 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 77 : < 0.04752007 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 78 : < 0.04752007 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 79 : < 0.04699171 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 80 : < 0.04699171 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 81 : < 0.04699171 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 82 : < 0.04475224 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 83 : < 0.04475224 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 84 : < 0.04475224 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 85 : < 0.03194127 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 86 : < 0.03194127 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 87 : < 0.03194127 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 88 : < 0.02446379 > ( 32.1367 ) 1.218388 0.1648005 0.3921812 89 : < 0.04054028 > ( 22.8408 ) 1.115788 0.1779293 0.337085 90 : < 0.04054028 > ( 22.8408 ) 1.115788 0.1779293 0.337085 91 : < 0.04054028 > ( 22.8408 ) 1.115788 0.1779293 0.337085 92 : < 0.04054028 > ( 22.8408 ) 1.115788 0.1779293 0.337085 93 : < 0.03934312 > ( 22.8408 ) 1.115788 0.1779293 0.337085 94 : < 0.03934312 > ( 22.8408 ) 1.115788 0.1779293 0.337085 95 : < 0.03934312 > ( 22.8408 ) 1.115788 0.1779293 0.337085 96 : < 0.03934312 > ( 22.8408 ) 1.115788 0.1779293 0.337085 97 : < 0.03934312 > ( 22.8408 ) 1.115788 0.1779293 0.337085 98 : < 0.03905163 > ( 22.8408 ) 1.115788 0.1779293 0.337085 99 : < 0.03354782 > ( 22.8408 ) 1.115788 0.1779293 0.337085 100 : < 0.03316904 > ( 22.8408 ) 1.115788 0.1779293 0.337085 101 : < 0.03316904 > ( 22.8408 ) 1.115788 0.1779293 0.337085 102 : < 0.03129767 > ( 22.8408 ) 1.115788 0.1779293 0.337085 103 : < 0.03129767 > ( 22.8408 ) 1.115788 0.1779293 0.337085 104 : < 0.02961744 > ( 22.8408 ) 1.115788 0.1779293 0.337085 105 : < 0.02961744 > ( 22.8408 ) 1.115788 0.1779293 0.337085 106 : < 0.02771651 > ( 22.8408 ) 1.115788 0.1779293 0.337085 107 : < 0.02771651 > ( 22.8408 ) 1.115788 0.1779293 0.337085 108 : < 0.02771651 > ( 22.8408 ) 1.115788 0.1779293 0.337085 109 : < 0.02771651 > ( 22.8408 ) 1.115788 0.1779293 0.337085 110 : < 0.02616082 > ( 22.8408 ) 1.115788 0.1779293 0.337085 111 : < 0.02616082 > ( 22.8408 ) 1.115788 0.1779293 0.337085 112 : < 0.02616082 > ( 22.8408 ) 1.115788 0.1779293 0.337085 113 : < 0.02616082 > ( 22.8408 ) 1.115788 0.1779293 0.337085 114 : < 0.02616082 > ( 22.8408 ) 1.115788 0.1779293 0.337085 115 : < 0.02616082 > ( 22.8408 ) 1.115788 0.1779293 0.337085 116 : < 0.02256137 > ( 22.8408 ) 1.115788 0.1779293 0.337085 117 : < 0.02250128 > ( 22.8408 ) 1.115788 0.1779293 0.337085 118 : < 0.02250128 > ( 22.8408 ) 1.115788 0.1779293 0.337085 119 : < 0.02212993 > ( 22.8408 ) 1.115788 0.1779293 0.337085 120 : < 0.01939336 > ( 22.0985 ) 0.9117826 0.2232727 0.2400501 121 : < 0.01935736 > ( 22.0985 ) 0.9117826 0.2232727 0.2400501 122 : < 0.01935736 > ( 22.0985 ) 0.9117826 0.2232727 0.2400501 123 : < 0.01935736 > ( 22.0985 ) 0.9117826 0.2232727 0.2400501 124 : < 0.01935736 > ( 22.0985 ) 0.9117826 0.2232727 0.2400501 125 : < 0.01935736 > ( 22.0985 ) 0.9117826 0.2232727 0.2400501 126 : < 0.01935736 > ( 22.0985 ) 0.9117826 0.2232727 0.2400501 127 : < 0.01935736 > ( 22.0985 ) 0.9117826 0.2232727 0.2400501 128 : < 0.01935736 > ( 22.0985 ) 0.9117826 0.2232727 0.2400501 129 : < 0.01753271 > ( 22.0985 ) 0.9117826 0.2232727 0.2400501 130 : < 0.02536802 > ( 17.56789 ) 1.115788 0.1820902 0.2877737 131 : < 0.02536802 > ( 17.56789 ) 1.115788 0.1820902 0.2877737 132 : < 0.02536802 > ( 17.56789 ) 1.115788 0.1820902 0.2877737 133 : < 0.02536802 > ( 17.56789 ) 1.115788 0.1820902 0.2877737 134 : < 0.02536802 > ( 17.56789 ) 1.115788 0.1820902 0.2877737 135 : < 0.02536802 > ( 17.56789 ) 1.115788 0.1820902 0.2877737 136 : < 0.01870163 > ( 17.56789 ) 1.115788 0.1820902 0.2877737 137 : < 0.01870163 > ( 17.56789 ) 1.115788 0.1820902 0.2877737 138 : < 0.01870163 > ( 17.56789 ) 1.115788 0.1820902 0.2877737 139 : < 0.01870163 > ( 17.56789 ) 1.115788 0.1820902 0.2877737 140 : < 0.01882865 > ( 17.49444 ) 1.115788 0.1820902 0.2886876 141 : < 0.01698791 > ( 17.49444 ) 1.115788 0.1820902 0.2886876 142 : < 0.01689915 > ( 17.49444 ) 1.115788 0.1820902 0.2886876 143 : < 0.01343627 > ( 17.49444 ) 1.115788 0.1820902 0.2886876 144 : < 0.01343627 > ( 17.49444 ) 1.115788 0.1820902 0.2886876 145 : < 0.01343627 > ( 17.49444 ) 1.115788 0.1820902 0.2886876 146 : < 0.01187015 > ( 17.49444 ) 1.115788 0.1820902 0.2886876 147 : < 0.01172604 > ( 17.49444 ) 1.115788 0.1820902 0.2886876 148 : < 0.01172604 > ( 17.49444 ) 1.115788 0.1820902 0.2886876 149 : < 0.01172604 > ( 17.49444 ) 1.115788 0.1820902 0.2886876 150 : < 0.01670415 > ( 14.61595 ) 0.9567954 0.2043113 0.2429337 151 : < 0.01670415 > ( 14.61595 ) 0.9567954 0.2043113 0.2429337 152 : < 0.01670415 > ( 14.61595 ) 0.9567954 0.2043113 0.2429337 153 : < 0.01670415 > ( 14.61595 ) 0.9567954 0.2043113 0.2429337 154 : < 0.01670415 > ( 14.61595 ) 0.9567954 0.2043113 0.2429337 155 : < 0.01670415 > ( 14.61595 ) 0.9567954 0.2043113 0.2429337 156 : < 0.01670415 > ( 14.61595 ) 0.9567954 0.2043113 0.2429337 157 : < 0.01670415 > ( 14.61595 ) 0.9567954 0.2043113 0.2429337 158 : < 0.01417443 > ( 14.61595 ) 0.9567954 0.2043113 0.2429337 159 : < 0.01417443 > ( 14.61595 ) 0.9567954 0.2043113 0.2429337 160 : < 0.01417443 > ( 14.61595 ) 0.9567954 0.2043113 0.2429337 161 : < 0.01417443 > ( 14.61595 ) 0.9567954 0.2043113 0.2429337 162 : < 0.01417443 > ( 14.61595 ) 0.9567954 0.2043113 0.2429337 163 : < 0.01417443 > ( 14.61595 ) 0.9567954 0.2043113 0.2429337 164 : < 0.01417443 > ( 14.61595 ) 0.9567954 0.2043113 0.2429337 165 : < 0.01417443 > ( 14.61595 ) 0.9567954 0.2043113 0.2429337 166 : < 0.01417443 > ( 14.61595 ) 0.9567954 0.2043113 0.2429337 167 : < 0.01240222 > ( 14.61595 ) 0.9567954 0.2043113 0.2429337 168 : < 0.01240222 > ( 14.61595 ) 0.9567954 0.2043113 0.2429337 169 : < 0.01240222 > ( 14.61595 ) 0.9567954 0.2043113 0.2429337 170 : < 0.01240222 > ( 14.61595 ) 0.9567954 0.2043113 0.2429337 171 : < 0.01240222 > ( 14.61595 ) 0.9567954 0.2043113 0.2429337 172 : < 0.01240222 > ( 14.61595 ) 0.9567954 0.2043113 0.2429337 173 : < 0.01240222 > ( 14.61595 ) 0.9567954 0.2043113 0.2429337 174 : < 0.01267251 > ( 14.45966 ) 0.9567954 0.2043113 0.3097119 175 : < 0.01267251 > ( 14.45966 ) 0.9567954 0.2043113 0.3097119 176 : < 0.01072199 > ( 14.45966 ) 0.9567954 0.2043113 0.3097119 177 : < 0.01460422 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 178 : < 0.01460422 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 179 : < 0.01460422 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 180 : < 0.01460422 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 181 : < 0.01460422 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 182 : < 0.01460422 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 183 : < 0.01460422 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 184 : < 0.01453987 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 185 : < 0.01450912 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 186 : < 0.01450912 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 187 : < 0.0141267 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 188 : < 0.0141267 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 189 : < 0.0141267 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 190 : < 0.0141267 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 191 : < 0.01283042 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 192 : < 0.01283042 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 193 : < 0.01283042 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 194 : < 0.01217961 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 195 : < 0.01217961 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 196 : < 0.01217961 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 197 : < 0.01217961 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 198 : < 0.01217961 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 199 : < 0.01133559 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 200 : < 0.01067607 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 201 : < 0.006848681 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 202 : < 0.006848681 > ( 12.21484 ) 0.9990202 0.1975343 0.2905052 203 : < 0.008647453 > ( 10.55097 ) 1.125506 0.1834127 0.2941167 204 : < 0.008647453 > ( 10.55097 ) 1.125506 0.1834127 0.2941167 205 : < 0.008647453 > ( 10.55097 ) 1.125506 0.1834127 0.2941167 206 : < 0.008647453 > ( 10.55097 ) 1.125506 0.1834127 0.2941167 207 : < 0.007927429 > ( 10.55097 ) 1.125506 0.1834127 0.2941167 208 : < 0.007927429 > ( 10.55097 ) 1.125506 0.1834127 0.2941167 209 : < 0.007769746 > ( 10.55097 ) 1.125506 0.1834127 0.2941167 210 : < 0.007769746 > ( 10.55097 ) 1.125506 0.1834127 0.2941167 211 : < 0.007769746 > ( 10.55097 ) 1.125506 0.1834127 0.2941167 212 : < 0.007769746 > ( 10.55097 ) 1.125506 0.1834127 0.2941167 213 : < 0.007769746 > ( 10.55097 ) 1.125506 0.1834127 0.2941167 214 : < 0.006777771 > ( 10.09105 ) 1.017428 0.197739 0.2636398 215 : < 0.006777771 > ( 10.09105 ) 1.017428 0.197739 0.2636398 216 : < 0.006722801 > ( 10.09105 ) 1.017428 0.197739 0.2636398 217 : < 0.006722801 > ( 10.09105 ) 1.017428 0.197739 0.2636398 218 : < 0.006722801 > ( 10.09105 ) 1.017428 0.197739 0.2636398 219 : < 0.006722801 > ( 10.09105 ) 1.017428 0.197739 0.2636398 220 : < 0.006722801 > ( 10.09105 ) 1.017428 0.197739 0.2636398 221 : < 0.006722801 > ( 10.09105 ) 1.017428 0.197739 0.2636398 222 : < 0.006722801 > ( 10.09105 ) 1.017428 0.197739 0.2636398 223 : < 0.006160303 > ( 10.09105 ) 1.017428 0.197739 0.2636398 224 : < 0.006160303 > ( 10.09105 ) 1.017428 0.197739 0.2636398 225 : < 0.006160303 > ( 10.09105 ) 1.017428 0.197739 0.2636398 226 : < 0.006160303 > ( 10.09105 ) 1.017428 0.197739 0.2636398 227 : < 0.005169555 > ( 10.09105 ) 1.017428 0.197739 0.2636398 228 : < 0.004996161 > ( 10.09105 ) 1.017428 0.197739 0.2636398 229 : < 0.004996161 > ( 10.09105 ) 1.017428 0.197739 0.2636398 230 : < 0.004996161 > ( 10.09105 ) 1.017428 0.197739 0.2636398 231 : < 0.004996161 > ( 10.09105 ) 1.017428 0.197739 0.2636398 232 : < 0.005084063 > ( 10.04022 ) 1.120079 0.184461 0.2915502 233 : < 0.003868612 > ( 10.04022 ) 1.120079 0.184461 0.2915502 234 : < 0.003226676 > ( 10.04022 ) 1.120079 0.184461 0.2915502 235 : < 0.003315432 > ( 9.988902 ) 1.017428 0.197739 0.2655706 236 : < 0.002828918 > ( 9.988902 ) 1.017428 0.197739 0.2655706 237 : < 0.002828918 > ( 9.988902 ) 1.017428 0.197739 0.2655706 238 : < 0.002828918 > ( 9.988902 ) 1.017428 0.197739 0.2655706 239 : < 0.002654422 > ( 9.988902 ) 1.017428 0.197739 0.2655706 240 : < 0.002646546 > ( 9.988902 ) 1.017428 0.197739 0.2655706 241 : < 0.002646546 > ( 9.988902 ) 1.017428 0.197739 0.2655706 242 : < 0.002551065 > ( 9.988902 ) 1.017428 0.197739 0.2655706 243 : < 0.002287334 > ( 9.988902 ) 1.017428 0.197739 0.2655706 244 : < 0.002094888 > ( 9.988902 ) 1.017428 0.197739 0.2655706 245 : < 0.001230538 > ( 9.988902 ) 1.017428 0.197739 0.2655706 246 : < 0.001230538 > ( 9.988902 ) 1.017428 0.197739 0.2655706 247 : < 0.001421687 > ( 9.728904 ) 0.9770017 0.2027812 0.2929567 248 : < 0.001421687 > ( 9.728904 ) 0.9770017 0.2027812 0.2929567 249 : < 0.001492989 > ( 9.687675 ) 0.9820805 0.202227 0.2857454 250 : < 0.001492989 > ( 9.687675 ) 0.9820805 0.202227 0.2857454 251 : < 0.001413056 > ( 9.687675 ) 0.9820805 0.202227 0.2857454 252 : < 0.001548094 > ( 9.545152 ) 0.9910307 0.2012494 0.2845529 253 : < 0.001896093 > ( 9.258389 ) 0.9876406 0.201579 0.2949691 254 : < 0.001490536 > ( 9.178351 ) 0.9928983 0.200902 0.2868676 255 : < 0.001368037 > ( 9.133262 ) 1.001906 0.1997992 0.2842197 256 : < 0.001287007 > ( 9.051984 ) 1.001626 0.1998094 0.2855302 257 : < 0.00110684 > ( 9.051984 ) 1.001626 0.1998094 0.2855302 258 : < 0.001367592 > ( 8.896893 ) 0.9951674 0.2004634 0.2842197 259 : < 0.0009684546 > ( 8.896893 ) 0.9951674 0.2004634 0.2842197 260 : < 0.0008341594 > ( 8.896893 ) 0.9951674 0.2004634 0.2842197 261 : < 0.0004344251 > ( 8.896893 ) 0.9951674 0.2004634 0.2842197 262 : < 0.0003967163 > ( 8.896893 ) 0.9951674 0.2004634 0.2842197 263 : < 0.0002943934 > ( 8.887488 ) 1.000065 0.1999005 0.284762 264 : < 0.0002382285 > ( 8.885436 ) 1.002874 0.1994934 0.2806912 265 : < 0.000173457 > ( 8.885436 ) 1.002874 0.1994934 0.2806912 266 : < 0.000160741 > ( 8.885436 ) 1.002874 0.1994934 0.2806912 267 : < 0.0001481587 > ( 8.885436 ) 1.002874 0.1994934 0.2806912 268 : < 0.0001168907 > ( 8.879155 ) 1.000572 0.199858 0.2865726 269 : < 8.345589e-05 > ( 8.879155 ) 1.000572 0.199858 0.2865726 270 : < 7.179479e-05 > ( 8.879155 ) 1.000572 0.199858 0.2865726 271 : < 0.0001339374 > ( 8.842449 ) 0.9964738 0.2002777 0.2837573 272 : < 0.0001033275 > ( 8.842449 ) 0.9964738 0.2002777 0.2837573 273 : < 0.0001410324 > ( 8.820646 ) 0.9991624 0.1999167 0.2813687 274 : < 0.0001350989 > ( 8.820646 ) 0.9991624 0.1999167 0.2813687 275 : < 0.0001255862 > ( 8.820646 ) 0.9991624 0.1999167 0.2813687 276 : < 0.0001182667 > ( 8.820646 ) 0.9991624 0.1999167 0.2813687 277 : < 0.000112048 > ( 8.820646 ) 0.9991624 0.1999167 0.2813687 278 : < 9.591002e-05 > ( 8.820646 ) 0.9991624 0.1999167 0.2813687 279 : < 8.075011e-05 > ( 8.820646 ) 0.9991624 0.1999167 0.2813687 280 : < 6.287093e-05 > ( 8.815232 ) 1.0002 0.1997529 0.2792734 281 : < 6.287093e-05 > ( 8.815232 ) 1.0002 0.1997529 0.2792734 282 : < 6.188651e-05 > ( 8.815232 ) 1.0002 0.1997529 0.2792734 283 : < 4.674904e-05 > ( 8.813412 ) 1.000361 0.1997252 0.2788947 284 : < 4.25444e-05 > ( 8.813412 ) 1.000361 0.1997252 0.2788947 285 : < 4.25444e-05 > ( 8.813412 ) 1.000361 0.1997252 0.2788947 286 : < 4.25444e-05 > ( 8.813412 ) 1.000361 0.1997252 0.2788947 287 : < 3.814525e-05 > ( 8.813412 ) 1.000361 0.1997252 0.2788947 288 : < 3.566861e-05 > ( 8.813412 ) 1.000361 0.1997252 0.2788947 289 : < 2.034132e-05 > ( 8.812989 ) 0.9984705 0.1999675 0.2799322 290 : < 1.82665e-05 > ( 8.809975 ) 0.9995279 0.1998427 0.2799926 291 : < 1.382387e-05 > ( 8.809975 ) 0.9995279 0.1998427 0.2799926 292 : < 1.339289e-05 > ( 8.809975 ) 0.9995279 0.1998427 0.2799926 293 : < 1.086682e-05 > ( 8.809975 ) 0.9995279 0.1998427 0.2799926 294 : < 1.142041e-05 > ( 8.809003 ) 0.9989595 0.1999131 0.2802209 295 : < 1.308622e-05 > ( 8.807544 ) 0.9987925 0.1998992 0.2785015 296 : < 1.308622e-05 > ( 8.807544 ) 0.9987925 0.1998992 0.2785015 297 : < 1.147405e-05 > ( 8.807544 ) 0.9987925 0.1998992 0.2785015 298 : < 1.128575e-05 > ( 8.807544 ) 0.9987925 0.1998992 0.2785015 299 : < 1.003239e-05 > ( 8.807544 ) 0.9987925 0.1998992 0.2785015 300 : < 9.416403e-06 > ( 8.807544 ) 0.9987925 0.1998992 0.2785015 301 : < 7.084855e-06 > ( 8.807544 ) 0.9987925 0.1998992 0.2785015 302 : < 6.339146e-06 > ( 8.807544 ) 0.9987925 0.1998992 0.2785015 303 : < 5.680015e-06 > ( 8.807544 ) 0.9987925 0.1998992 0.2785015 304 : < 5.680015e-06 > ( 8.807544 ) 0.9987925 0.1998992 0.2785015 305 : < 5.551753e-06 > ( 8.807544 ) 0.9987925 0.1998992 0.2785015 306 : < 4.264251e-06 > ( 8.807544 ) 0.9987925 0.1998992 0.2785015 307 : < 4.264251e-06 > ( 8.807544 ) 0.9987925 0.1998992 0.2785015 308 : < 4.264251e-06 > ( 8.807544 ) 0.9987925 0.1998992 0.2785015 309 : < 3.310676e-06 > ( 8.807544 ) 0.9987925 0.1998992 0.2785015 310 : < 2.90051e-06 > ( 8.807544 ) 0.9987925 0.1998992 0.2785015 311 : < 2.868159e-06 > ( 8.807544 ) 0.9987925 0.1998992 0.2785015 312 : < 2.578554e-06 > ( 8.807496 ) 0.9991107 0.1998889 0.2799183 313 : < 3.251148e-06 > ( 8.806875 ) 0.9988974 0.1999144 0.2799686 314 : < 3.178327e-06 > ( 8.806875 ) 0.9988974 0.1999144 0.2799686 315 : < 2.097213e-06 > ( 8.806875 ) 0.9988974 0.1999144 0.2799686 316 : < 1.48446e-06 > ( 8.806875 ) 0.9988974 0.1999144 0.2799686 317 : < 1.155885e-06 > ( 8.806875 ) 0.9988974 0.1999144 0.2799686 318 : < 1.116642e-06 > ( 8.806875 ) 0.9988974 0.1999144 0.2799686 319 : < 1.116642e-06 > ( 8.806875 ) 0.9988974 0.1999144 0.2799686 320 : < 1.071415e-06 > ( 8.806875 ) 0.9988974 0.1999144 0.2799686 321 : < 1.071415e-06 > ( 8.806875 ) 0.9988974 0.1999144 0.2799686 322 : < 9.53625e-07 > ( 8.806875 ) 0.9988974 0.1999144 0.2799686 323 : < 9.53625e-07 > ( 8.806875 ) 0.9988974 0.1999144 0.2799686 324 : < 7.139284e-07 > ( 8.806875 ) 0.9988974 0.1999144 0.2799686 325 : < 5.401451e-07 > ( 8.806875 ) 0.9988974 0.1999144 0.2799686 326 : < 7.729068e-07 > ( 8.806741 ) 0.9988289 0.1999198 0.2798315 327 : < 1.01155e-06 > ( 8.806595 ) 0.9987117 0.1999334 0.2798342 328 : < 9.327248e-07 > ( 8.806595 ) 0.9987117 0.1999334 0.2798342 329 : < 8.305895e-07 > ( 8.806571 ) 0.99848 0.1999564 0.2796244 330 : < 8.276247e-07 > ( 8.806553 ) 0.9987529 0.1999192 0.2793319 331 : < 3.658393e-07 > ( 8.806529 ) 0.9984135 0.1999515 0.2789616 332 : < 7.33075e-07 > ( 8.806257 ) 0.9986107 0.1999354 0.2793235 333 : < 5.856533e-07 > ( 8.806257 ) 0.9986107 0.1999354 0.2793235 334 : < 5.856533e-07 > ( 8.806257 ) 0.9986107 0.1999354 0.2793235 335 : < 5.70618e-07 > ( 8.806257 ) 0.9986107 0.1999354 0.2793235 336 : < 5.214147e-07 > ( 8.806257 ) 0.9986107 0.1999354 0.2793235 337 : < 5.214147e-07 > ( 8.806257 ) 0.9986107 0.1999354 0.2793235 338 : < 5.097635e-07 > ( 8.806257 ) 0.9986107 0.1999354 0.2793235 339 : < 5.071348e-07 > ( 8.806257 ) 0.9986107 0.1999354 0.2793235 340 : < 3.215318e-07 > ( 8.806257 ) 0.9986107 0.1999354 0.2793235 341 : < 2.525058e-07 > ( 8.806257 ) 0.9986107 0.1999354 0.2793235 342 : < 2.478916e-07 > ( 8.806242 ) 0.9985037 0.1999447 0.2791625 343 : < 3.36089e-07 > ( 8.806168 ) 0.9987386 0.1999197 0.2792836 344 : < 3.232179e-07 > ( 8.806168 ) 0.9987386 0.1999197 0.2792836 345 : < 2.29253e-07 > ( 8.806168 ) 0.9987386 0.1999197 0.2792836 346 : < 2.242143e-07 > ( 8.806168 ) 0.9987386 0.1999197 0.2792836 347 : < 1.900574e-07 > ( 8.806168 ) 0.9987386 0.1999197 0.2792836 348 : < 1.786672e-07 > ( 8.806168 ) 0.9987386 0.1999197 0.2792836 349 : < 1.535724e-07 > ( 8.806168 ) 0.9987386 0.1999197 0.2792836 350 : < 1.05253e-07 > ( 8.806168 ) 0.9987386 0.1999197 0.2792836 351 : < 8.99405e-08 > ( 8.806168 ) 0.9987386 0.1999197 0.2792836 Time user system elapsed Time 1.02 0.01 1.03 > showProc.time() Time (user system elapsed): 2.4 0.01 2.42 > > plot(y ~ x, d.exp30, main = "Data = d.exp30") > cTr <- adjustcolor("red4", 0.5) > cLS <- adjustcolor("blue2", 0.5) > cE <- curve(Expo(x, a=1, b=0.2), 0, 10, n=1+2^9, col=cTr, lwd=2, lty=2, add=TRUE) > lines(d.exp30$x, fitted(Cfit), col=cLS, lwd=3) > ll <- length(m1 <- sapply(ls.str(patt="^Rfit"), get, simplify=FALSE)) > .tmp <- lapply(m1, function(.) lines(d.exp30$x, fitted(.))) > legend("topleft", c("true", "LS", names(m1)), + lwd=c(2,3, rep(1,ll)), lty=c(2,1, rep(1,ll)), + col=c(cTr,cLS, rep(par("fg"),ll)), bty="n", inset=.01) > showProc.time() Time (user system elapsed): 0 0 0 > > ## 40% outliers present {use different data name: seen in print() > d.exp40out <- within(d.exp30, y[15:27] <- y[15:27] + 100) > op <- options(digits=12) > Cfit.40out <- update(Cfit, data = d.exp40out, trace=TRUE, + control = nls.control(tol = Cfit$control$tol)) 130520.368888 (1.09e+00): par = (1 0.2) 80939.4887213 (6.97e-01): par = (3.64667878025 0.0101134659897) 62863.1807960 (9.98e-02): par = (3.17512533113 0.139117100179) 62293.4269296 (2.59e-02): par = (3.33096092124 0.127875736076) 62277.8108522 (1.65e-02): par = (3.27404212048 0.136208504237) 62272.2383633 (1.11e-02): par = (3.30909973778 0.130726117097) 62269.6374531 (7.32e-03): par = (3.28557788203 0.134382004845) 62268.5303739 (4.89e-03): par = (3.30112850767 0.131955807605) 62268.0284036 (3.24e-03): par = (3.29074492291 0.133571694035) 62267.8100195 (2.16e-03): par = (3.2976336361 0.13249786192) 62267.7123190 (1.44e-03): par = (3.29304349577 0.133212576944) 62267.6694102 (9.56e-04): par = (3.29609320105 0.132737361813) 62267.6503306 (6.35e-04): par = (3.29406308865 0.133053542914) 62267.6419171 (4.23e-04): par = (3.29541273769 0.132843271135) 62267.6381871 (2.81e-04): par = (3.29451465226 0.132983159173) 62267.6365387 (1.87e-04): par = (3.29511194287 0.132890110001) 62267.6358086 (1.24e-04): par = (3.29471456834 0.132952009247) 62267.6354859 (8.28e-05): par = (3.2949789053 0.132910830735) 62267.6353430 (5.51e-05): par = (3.29480300795 0.132938230727) 62267.6352797 (3.67e-05): par = (3.29492005295 0.132919998024) 62267.6352517 (2.44e-05): par = (3.29484215578 0.132932132464) 62267.6352393 (1.62e-05): par = (3.29489399679 0.132924056475) 62267.6352338 (1.08e-05): par = (3.29485953776 0.132929424915) 62267.6352314 (7.22e-06): par = (3.294882448 0.132925855577) 62267.6352303 (4.83e-06): par = (3.29486712277 0.132928243013) 62267.6352298 (3.20e-06): par = (3.2948773806 0.132926644989) 62267.6352296 (2.11e-06): par = (3.29487057599 0.132927704924) 62267.6352295 (1.39e-06): par = (3.29487506278 0.132927006299) 62267.6352295 (9.00e-07): par = (3.29487209824 0.132927467432) 62267.6352294 (5.95e-07): par = (3.29487401301 0.132927169617) 62267.6352294 (4.08e-07): par = (3.29487275137 0.132927366302) 62267.6352294 (2.60e-07): par = (3.29487361677 0.132927231339) 62267.6352294 (1.81e-07): par = (3.29487306512 0.132927317416) 62267.6352294 (1.25e-07): par = (3.29487344896 0.132927257691) 62267.6352294 (5.66e-08): par = (3.2948731801 0.132927299145) > if(FALSE) ## this fails for "bad" non-R BLAS/LAPACK + Cfit.no.out <- update(Cfit.40out, subset = -(15:27)) > ## help giving it a good start *and* raise tolerance (from 8e-8): > ## still fails for all three of {ATLAS, MKL, OpenBLAS} with > ## Error in nls(formula = y ~ Expo(x, a, b), data = d.exp.Hlev, start = c(a = 1, : > ## step factor 0.000488281 reduced below 'minFactor' of 0.000976562 > Cfit.no.out <- + tryCatch(error = function(e) e, + update(Cfit.40out, subset = -(15:27), start = c(a = 1, b = 0.2), trace=TRUE, + control = nls.control(maxiter = 1000, tol = 5e-7, printEval=TRUE)) + ) 11.1131221181 (2.38e-01): par = (1 0.2) It. 1, fac= 1, eval (no.,total): ( 1, 1): new dev = 10.5154 10.5153984590 (2.73e-03): par = (0.983651552028 0.199175843153) It. 2, fac= 1, eval (no.,total): ( 1, 2): new dev = 10.5153 10.5153201605 (1.22e-06): par = (0.983507348145 0.199160280981) It. 3, fac= 1, eval (no.,total): ( 1, 3): new dev = 10.5153 10.5153201604 (2.26e-09): par = (0.983507009529 0.19916032082) > Cfit.no..ok <- !inherits(Cfit.no.out, "error") > options(op) > if(doExtras) { + Rf.out.MM.S.bisquare <- update(Rfit.MM.S.bisquare, data=d.exp40out) + Rf.out.MM.S.lqq <- update(Rf.out.MM.S.bisquare, psi = "lqq") + Rf.out.MM.S.optimal <- update(Rf.out.MM.S.bisquare, psi = "optimal") + Rf.out.MM.S.hampel <- update(Rf.out.MM.S.bisquare, psi = "hampel") + showProc.time() + } > Rf.out.MM.lts.bisquare <- update(Rfit.MM.S.bisquare, data=d.exp40out, init= "lts") > Rf.out.MM.lts.lqq <- update(Rf.out.MM.lts.bisquare, psi= "lqq") #----------- > Rf.out.MM.lts.optimal <- update(Rf.out.MM.lts.bisquare, psi= "optimal") > Rf.out.MM.lts.hampel <- update(Rf.out.MM.lts.bisquare, psi= "hampel") > showProc.time() Time (user system elapsed): 0.36 0 0.36 > > Rf.out.tau.bisquare <- update(Rfit.tau.bisquare, data=d.exp40out) > Rf.out.tau.optimal <- update(Rfit.tau.bisquare, data=d.exp40out, psi = "optimal") > Rf.out.CM <- update(Rfit.CM, data=d.exp40out) > Rf.out.mtl <- update(Rfit.mtl, data=d.exp40out) 1 : < 0.003562711 > ( 42.35703 ) 1 0.2 1 2 : < 0.002683904 > ( 42.35703 ) 1 0.2 1 3 : < 0.002299342 > ( 42.35703 ) 1 0.2 1 4 : < 0.00208893 > ( 42.35703 ) 1 0.2 1 5 : < 0.002040813 > ( 42.35703 ) 1 0.2 1 6 : < 0.00197901 > ( 42.35703 ) 1 0.2 1 7 : < 0.001938679 > ( 42.35703 ) 1 0.2 1 8 : < 0.001875206 > ( 42.35703 ) 1 0.2 1 9 : < 0.001718907 > ( 42.35703 ) 1 0.2 1 10 : < 0.00143452 > ( 42.35703 ) 1 0.2 1 11 : < 0.001293953 > ( 42.35703 ) 1 0.2 1 12 : < 0.001167517 > ( 42.35703 ) 1 0.2 1 13 : < 0.001147897 > ( 42.35703 ) 1 0.2 1 14 : < 0.001022006 > ( 42.35703 ) 1 0.2 1 15 : < 0.001004856 > ( 42.35703 ) 1 0.2 1 16 : < 0.001004856 > ( 42.35703 ) 1 0.2 1 17 : < 0.001004856 > ( 42.35703 ) 1 0.2 1 18 : < 0.0009918255 > ( 42.35703 ) 1 0.2 1 19 : < 0.0009906996 > ( 42.35703 ) 1 0.2 1 20 : < 0.0008356166 > ( 42.35703 ) 1 0.2 1 21 : < 0.0007500565 > ( 42.35703 ) 1 0.2 1 22 : < 0.0007500565 > ( 42.35703 ) 1 0.2 1 23 : < 0.0007138254 > ( 42.35703 ) 1 0.2 1 24 : < 0.0006892772 > ( 42.35703 ) 1 0.2 1 25 : < 0.0006892772 > ( 42.35703 ) 1 0.2 1 26 : < 0.0006180502 > ( 42.35703 ) 1 0.2 1 27 : < 0.0004472576 > ( 42.35703 ) 1 0.2 1 28 : < 0.0004107399 > ( 42.35703 ) 1 0.2 1 29 : < 0.0004107399 > ( 42.35703 ) 1 0.2 1 30 : < 0.0004107399 > ( 42.35703 ) 1 0.2 1 31 : < 0.0003780804 > ( 42.35703 ) 1 0.2 1 32 : < 0.0003780804 > ( 42.35703 ) 1 0.2 1 33 : < 0.0003641482 > ( 42.35703 ) 1 0.2 1 34 : < 0.0003552235 > ( 42.35703 ) 1 0.2 1 35 : < 0.0003552235 > ( 42.35703 ) 1 0.2 1 36 : < 0.0003173703 > ( 42.35703 ) 1 0.2 1 37 : < 0.0003103385 > ( 42.35703 ) 1 0.2 1 38 : < 0.0003023294 > ( 41.2185 ) 1 0.2 0.8749227 39 : < 0.0003023294 > ( 41.2185 ) 1 0.2 0.8749227 40 : < 0.000285123 > ( 41.2185 ) 1 0.2 0.8749227 41 : < 0.000285123 > ( 41.2185 ) 1 0.2 0.8749227 42 : < 0.0002755209 > ( 41.2185 ) 1 0.2 0.8749227 43 : < 0.0002579781 > ( 41.2185 ) 1 0.2 0.8749227 44 : < 0.000248941 > ( 41.2185 ) 1 0.2 0.8749227 45 : < 0.0002087005 > ( 41.2185 ) 1 0.2 0.8749227 46 : < 0.0002087005 > ( 41.2185 ) 1 0.2 0.8749227 47 : < 0.0001798612 > ( 41.2185 ) 1 0.2 0.8749227 48 : < 0.0001737002 > ( 40.44839 ) 0.9757931 0.2 0.8749227 49 : < 0.0001631675 > ( 40.44839 ) 0.9757931 0.2 0.8749227 50 : < 0.0001546442 > ( 40.44839 ) 0.9757931 0.2 0.8749227 51 : < 0.0001346317 > ( 40.44839 ) 0.9757931 0.2 0.8749227 52 : < 0.0001292555 > ( 40.44839 ) 0.9757931 0.2 0.8749227 53 : < 0.0001243123 > ( 40.44839 ) 0.9757931 0.2 0.8749227 54 : < 0.000103482 > ( 40.44839 ) 0.9757931 0.2 0.8749227 55 : < 9.590413e-05 > ( 40.44839 ) 0.9757931 0.2 0.8749227 56 : < 8.887343e-05 > ( 40.44839 ) 0.9757931 0.2 0.8749227 57 : < 9.282321e-05 > ( 40.13592 ) 0.9757931 0.2 0.8163305 58 : < 9.282321e-05 > ( 40.13592 ) 0.9757931 0.2 0.8163305 59 : < 8.444278e-05 > ( 40.13592 ) 0.9757931 0.2 0.8163305 60 : < 5.399421e-05 > ( 40.13592 ) 0.9757931 0.2 0.8163305 61 : < 4.794108e-05 > ( 40.13592 ) 0.9757931 0.2 0.8163305 62 : < 4.138473e-05 > ( 40.13592 ) 0.9757931 0.2 0.8163305 63 : < 3.768665e-05 > ( 40.13592 ) 0.9757931 0.2 0.8163305 64 : < 3.768665e-05 > ( 40.13592 ) 0.9757931 0.2 0.8163305 65 : < 3.668443e-05 > ( 40.13592 ) 0.9757931 0.2 0.8163305 66 : < 3.31901e-05 > ( 40.13592 ) 0.9757931 0.2 0.8163305 67 : < 2.896316e-05 > ( 40.13592 ) 0.9757931 0.2 0.8163305 68 : < 2.204796e-05 > ( 40.13592 ) 0.9757931 0.2 0.8163305 69 : < 1.843594e-05 > ( 40.13592 ) 0.9757931 0.2 0.8163305 70 : < 1.635856e-05 > ( 40.13592 ) 0.9757931 0.2 0.8163305 71 : < 1.635856e-05 > ( 40.13592 ) 0.9757931 0.2 0.8163305 72 : < 1.415656e-05 > ( 40.13379 ) 0.9917102 0.1986594 0.8058799 73 : < 1.225375e-05 > ( 40.13379 ) 0.9917102 0.1986594 0.8058799 74 : < 9.543117e-06 > ( 40.13379 ) 0.9917102 0.1986594 0.8058799 75 : < 8.01748e-06 > ( 40.13379 ) 0.9917102 0.1986594 0.8058799 76 : < 4.502424e-06 > ( 40.10744 ) 0.9844145 0.1988434 0.8082025 77 : < 4.218635e-06 > ( 40.10744 ) 0.9844145 0.1988434 0.8082025 78 : < 1.823543e-06 > ( 40.10744 ) 0.9844145 0.1988434 0.8082025 79 : < 1.764368e-06 > ( 40.10744 ) 0.9844145 0.1988434 0.8082025 80 : < 1.395215e-06 > ( 40.10398 ) 0.9783252 0.1998861 0.7686048 81 : < 1.129953e-06 > ( 40.10398 ) 0.9783252 0.1998861 0.7686048 82 : < 1.375728e-06 > ( 40.08288 ) 0.9813502 0.1994292 0.7954821 83 : < 1.375728e-06 > ( 40.08288 ) 0.9813502 0.1994292 0.7954821 84 : < 1.375728e-06 > ( 40.08288 ) 0.9813502 0.1994292 0.7954821 85 : < 1.11625e-06 > ( 40.08288 ) 0.9813502 0.1994292 0.7954821 86 : < 7.473127e-07 > ( 40.08288 ) 0.9813502 0.1994292 0.7954821 87 : < 3.521701e-07 > ( 40.08288 ) 0.9813502 0.1994292 0.7954821 88 : < 3.104966e-07 > ( 40.08288 ) 0.9813502 0.1994292 0.7954821 89 : < 3.104966e-07 > ( 40.08288 ) 0.9813502 0.1994292 0.7954821 90 : < 1.929642e-07 > ( 40.08288 ) 0.9813502 0.1994292 0.7954821 91 : < 1.319867e-07 > ( 40.08288 ) 0.9813502 0.1994292 0.7954821 92 : < 1.707788e-07 > ( 40.07981 ) 0.9816332 0.199275 0.7891589 93 : < 1.230306e-07 > ( 40.07981 ) 0.9816332 0.199275 0.7891589 94 : < 1.108068e-07 > ( 40.07981 ) 0.9816332 0.199275 0.7891589 95 : < 1.108068e-07 > ( 40.07981 ) 0.9816332 0.199275 0.7891589 96 : < 8.860688e-08 > ( 40.07838 ) 0.9819907 0.1993515 0.7891589 > showProc.time() Time (user system elapsed): 1.52 0.02 1.54 > > plot(y ~ x, d.exp40out, main = "Data = d.exp40out") > cE <- curve(Expo(x, a=1, b=0.2), 0, 10, n=1+2^9, col=cTr, lwd=2, lty=2, add=TRUE) > ll <- length(m1 <- sapply(ls.str(patt="^Rf.out"), get, simplify=FALSE)) > .tmp <- lapply(m1, function(.) lines(d.exp40out$x, fitted(.))) > xx <- local({p <- par("usr"); seq(p[1],p[2], length.out=256)}) > if(Cfit.no..ok) + lines(xx, predict(Cfit.no.out, list(x=xx)), col=cLS, lwd=3) > lines(xx, predict(Cfit.40out , list(x=xx)), col=cLS, lty=2) > legend("topleft", c("true", "LS [w/o outl]", "LS", names(m1)), + lwd=c(2,3, rep(1,1+ll)), lty=c(2,1,2, rep(1,ll)), + col=c(cTr,cLS,cLS, rep(par("fg"),ll)), bty="n", inset=.01) > showProc.time() Time (user system elapsed): 0.01 0 0.01 > > ## presence of high leverage point outliers > d.exp.Hlev <- within(d.exp40out, { + x[28:30] <- x[28:30] + 10 ## shift 10 + y <- Expo(x, 1, 0.2) + err + y[28:30] <- y[28:30] + 500 + }) > > op <- options(digits=12) > Cfit.Hlev <- + tryCatch(error = function(e) e, + update(Cfit.40out, data = d.exp.Hlev, start = c(a = 1, b = 0.2), trace=TRUE, + control = nls.control(maxiter = 100, tol = 5e-7, printEval=TRUE)) + ) 748655.257776 (7.07e+00): par = (1 0.2) It. 1, fac= 1, eval (no.,total): ( 1, 1): new dev = 1.05162e+08 It. 1, fac= 0.5, eval (no.,total): ( 2, 2): new dev = 266799 266799.286509 (4.73e+00): par = (0.638594559239 0.321213457043) It. 2, fac= 1, eval (no.,total): ( 1, 3): new dev = 29382 29382.0352946 (9.18e-01): par = (2.32791139589 0.217182750239) It. 3, fac= 1, eval (no.,total): ( 1, 4): new dev = 17419.3 17419.2743066 (3.47e-02): par = (1.55875387852 0.254039889824) It. 4, fac= 1, eval (no.,total): ( 1, 5): new dev = 17398.6 17398.6417634 (3.03e-03): par = (1.54249909618 0.25467120046) It. 5, fac= 1, eval (no.,total): ( 1, 6): new dev = 17398.6 17398.5581694 (1.45e-03): par = (1.54815322201 0.254375395115) It. 6, fac= 1, eval (no.,total): ( 1, 7): new dev = 17398.5 17398.5389073 (6.89e-04): par = (1.54545160428 0.254516507295) It. 7, fac= 1, eval (no.,total): ( 1, 8): new dev = 17398.5 17398.5345908 (3.28e-04): par = (1.54673528957 0.254449455066) It. 8, fac= 1, eval (no.,total): ( 1, 9): new dev = 17398.5 17398.5336075 (1.56e-04): par = (1.54612393961 0.25448138788) It. 9, fac= 1, eval (no.,total): ( 1, 10): new dev = 17398.5 17398.5333856 (7.42e-05): par = (1.54641462827 0.254466204114) It. 10, fac= 1, eval (no.,total): ( 1, 11): new dev = 17398.5 17398.5333352 (3.52e-05): par = (1.54627646809 0.25447342072) It. 11, fac= 1, eval (no.,total): ( 1, 12): new dev = 17398.5 17398.5333239 (1.67e-05): par = (1.54634206606 0.254469994315) It. 12, fac= 1, eval (no.,total): ( 1, 13): new dev = 17398.5 17398.5333213 (7.82e-06): par = (1.54631098239 0.254471617924) It. 13, fac= 1, eval (no.,total): ( 1, 14): new dev = 17398.5 17398.5333208 (3.72e-06): par = (1.54632554926 0.254470857022) It. 14, fac= 1, eval (no.,total): ( 1, 15): new dev = 17398.5 17398.5333206 (1.73e-06): par = (1.54631862747 0.254471218589) It. 15, fac= 1, eval (no.,total): ( 1, 16): new dev = 17398.5 17398.5333206 (6.81e-07): par = (1.54632185654 0.254471049916) It. 16, fac= 1, eval (no.,total): ( 1, 17): new dev = 17398.5 17398.5333206 (2.64e-07): par = (1.54632058784 0.254471116231) > if(Cfit.Hlev..ok <- !inherits(Cfit.Hlev, "error")) { + Cfit.no.Hlev <- update(Cfit.Hlev, subset = -(28:30)) + } else { ## substitute -- better? + Cfit.no.Hlev <- update(Cfit, subset = -(28:30)) + } 18.7899627711 (3.66e-01): par = (1 0.2) It. 1, fac= 1, eval (no.,total): ( 1, 1): new dev = 16.5647 16.5646938715 (7.02e-03): par = (0.920095152037 0.216822360981) It. 2, fac= 1, eval (no.,total): ( 1, 2): new dev = 16.5639 16.5638786808 (3.36e-06): par = (0.919287416513 0.216838534591) It. 3, fac= 1, eval (no.,total): ( 1, 3): new dev = 16.5639 16.5638786806 (2.58e-08): par = (0.919286403768 0.216838659367) > showProc.time() Time (user system elapsed): 0 0 0 > options(op) > > if(doExtras) { + Rf.Hlev.MM.S.bisquare <- update(Rfit.MM.S.bisquare, data = d.exp.Hlev) + Rf.Hlev.MM.S.lqq <- update(Rf.Hlev.MM.S.bisquare, psi = "lqq") + Rf.Hlev.MM.S.optimal <- update(Rf.Hlev.MM.S.bisquare, psi = "optimal") + Rf.Hlev.MM.S.hampel <- update(Rf.Hlev.MM.S.bisquare, psi = "hampel") + showProc.time() + } > Rf.Hlev.MM.lts.bisquare <- update(Rfit.MM.S.bisquare, data = d.exp.Hlev, init="lts") > Rf.Hlev.MM.lts.lqq <- update(Rf.Hlev.MM.lts.bisquare, psi= "lqq") > Rf.Hlev.MM.lts.optimal <- update(Rf.Hlev.MM.lts.bisquare, psi="optimal") > Rf.Hlev.MM.lts.hampel <- update(Rf.Hlev.MM.lts.bisquare, psi= "hampel") > showProc.time() Time (user system elapsed): 0.36 0 0.36 > > Rf.Hlev.tau.bisquare <- update(Rfit.tau.bisquare, data = d.exp.Hlev) > Rf.Hlev.tau.optimal <- update(Rf.Hlev.tau.bisquare, psi = "optimal") > Rf.Hlev.CM <- update(Rfit.CM, data = d.exp.Hlev) > Rf.Hlev.mtl <- update(Rfit.mtl, data = d.exp.Hlev) 1 : < 9.932221e-05 > ( 68.41264 ) 1 0.2 1 2 : < 6.274736e-05 > ( 68.41264 ) 1 0.2 1 3 : < 6.274736e-05 > ( 68.41264 ) 1 0.2 1 4 : < 5.392995e-05 > ( 68.41264 ) 1 0.2 1 5 : < 5.088127e-05 > ( 68.41264 ) 1 0.2 1 6 : < 4.908561e-05 > ( 68.41264 ) 1 0.2 1 7 : < 4.605234e-05 > ( 68.41264 ) 1 0.2 1 8 : < 4.47368e-05 > ( 68.41264 ) 1 0.2 1 9 : < 4.415354e-05 > ( 68.41264 ) 1 0.2 1 10 : < 4.368908e-05 > ( 68.41264 ) 1 0.2 1 11 : < 4.110801e-05 > ( 68.41264 ) 1 0.2 1 12 : < 4.069338e-05 > ( 68.41264 ) 1 0.2 1 13 : < 4.031166e-05 > ( 68.41264 ) 1 0.2 1 14 : < 3.892498e-05 > ( 65.18417 ) 1.189142 0.1262824 1 15 : < 3.832165e-05 > ( 65.18417 ) 1.189142 0.1262824 1 16 : < 3.646413e-05 > ( 65.18417 ) 1.189142 0.1262824 1 17 : < 3.346548e-05 > ( 65.18417 ) 1.189142 0.1262824 1 18 : < 3.202559e-05 > ( 65.18417 ) 1.189142 0.1262824 1 19 : < 3.202559e-05 > ( 65.18417 ) 1.189142 0.1262824 1 20 : < 3.202559e-05 > ( 65.18417 ) 1.189142 0.1262824 1 21 : < 3.144177e-05 > ( 65.18417 ) 1.189142 0.1262824 1 22 : < 3.144177e-05 > ( 65.18417 ) 1.189142 0.1262824 1 23 : < 2.902388e-05 > ( 57.94796 ) 1.71134 0.017981 1.041539 24 : < 2.757959e-05 > ( 57.94796 ) 1.71134 0.017981 1.041539 25 : < 3.265521e-05 > ( 52.00035 ) 1.597033 0.04463572 0.8763893 26 : < 3.230138e-05 > ( 52.00035 ) 1.597033 0.04463572 0.8763893 27 : < 3.227088e-05 > ( 52.00035 ) 1.597033 0.04463572 0.8763893 28 : < 3.165079e-05 > ( 52.00035 ) 1.597033 0.04463572 0.8763893 29 : < 2.850662e-05 > ( 52.00035 ) 1.597033 0.04463572 0.8763893 30 : < 2.327611e-05 > ( 52.00035 ) 1.597033 0.04463572 0.8763893 31 : < 2.327611e-05 > ( 52.00035 ) 1.597033 0.04463572 0.8763893 32 : < 3.5326e-05 > ( 39.43062 ) 0.4774236 0.3077108 0.5941036 33 : < 3.385232e-05 > ( 37.46558 ) 1.366756 0.1018557 0.6666641 34 : < 3.385232e-05 > ( 37.46558 ) 1.366756 0.1018557 0.6666641 35 : < 3.299797e-05 > ( 37.46558 ) 1.366756 0.1018557 0.6666641 36 : < 3.299797e-05 > ( 37.46558 ) 1.366756 0.1018557 0.6666641 37 : < 3.290589e-05 > ( 37.46558 ) 1.366756 0.1018557 0.6666641 38 : < 3.701338e-05 > ( 30.15308 ) 1.201137 0.1540887 0.4810447 39 : < 3.701338e-05 > ( 30.15308 ) 1.201137 0.1540887 0.4810447 40 : < 3.645266e-05 > ( 30.15308 ) 1.201137 0.1540887 0.4810447 41 : < 3.645266e-05 > ( 30.15308 ) 1.201137 0.1540887 0.4810447 42 : < 3.645266e-05 > ( 30.15308 ) 1.201137 0.1540887 0.4810447 43 : < 3.645266e-05 > ( 30.15308 ) 1.201137 0.1540887 0.4810447 44 : < 3.645266e-05 > ( 30.15308 ) 1.201137 0.1540887 0.4810447 45 : < 3.087764e-05 > ( 30.15308 ) 1.201137 0.1540887 0.4810447 46 : < 3.087764e-05 > ( 30.15308 ) 1.201137 0.1540887 0.4810447 47 : < 3.087764e-05 > ( 30.15308 ) 1.201137 0.1540887 0.4810447 48 : < 3.087764e-05 > ( 30.15308 ) 1.201137 0.1540887 0.4810447 49 : < 3.087764e-05 > ( 30.15308 ) 1.201137 0.1540887 0.4810447 50 : < 3.087764e-05 > ( 30.15308 ) 1.201137 0.1540887 0.4810447 51 : < 3.087764e-05 > ( 30.15308 ) 1.201137 0.1540887 0.4810447 52 : < 2.657404e-05 > ( 30.15308 ) 1.201137 0.1540887 0.4810447 53 : < 3.096098e-05 > ( 25.63392 ) 1.122918 0.1764149 0.3855222 54 : < 3.096098e-05 > ( 25.63392 ) 1.122918 0.1764149 0.3855222 55 : < 3.096098e-05 > ( 25.63392 ) 1.122918 0.1764149 0.3855222 56 : < 3.096098e-05 > ( 25.63392 ) 1.122918 0.1764149 0.3855222 57 : < 3.096098e-05 > ( 25.63392 ) 1.122918 0.1764149 0.3855222 58 : < 3.096098e-05 > ( 25.63392 ) 1.122918 0.1764149 0.3855222 59 : < 3.096098e-05 > ( 25.63392 ) 1.122918 0.1764149 0.3855222 60 : < 3.096098e-05 > ( 25.63392 ) 1.122918 0.1764149 0.3855222 61 : < 2.744674e-05 > ( 25.63392 ) 1.122918 0.1764149 0.3855222 62 : < 2.744674e-05 > ( 25.63392 ) 1.122918 0.1764149 0.3855222 63 : < 2.744674e-05 > ( 25.63392 ) 1.122918 0.1764149 0.3855222 64 : < 2.394768e-05 > ( 25.63392 ) 1.122918 0.1764149 0.3855222 65 : < 2.394768e-05 > ( 25.63392 ) 1.122918 0.1764149 0.3855222 66 : < 2.394768e-05 > ( 25.63392 ) 1.122918 0.1764149 0.3855222 67 : < 2.785765e-05 > ( 21.60612 ) 1.049891 0.1896514 0.345301 68 : < 2.696136e-05 > ( 21.60612 ) 1.049891 0.1896514 0.345301 69 : < 2.696136e-05 > ( 21.60612 ) 1.049891 0.1896514 0.345301 70 : < 2.696136e-05 > ( 21.60612 ) 1.049891 0.1896514 0.345301 71 : < 2.236592e-05 > ( 21.60612 ) 1.049891 0.1896514 0.345301 72 : < 2.236592e-05 > ( 21.60612 ) 1.049891 0.1896514 0.345301 73 : < 2.236592e-05 > ( 21.60612 ) 1.049891 0.1896514 0.345301 74 : < 1.399105e-05 > ( 21.60612 ) 1.049891 0.1896514 0.345301 75 : < 1.399105e-05 > ( 21.60612 ) 1.049891 0.1896514 0.345301 76 : < 1.399105e-05 > ( 21.60612 ) 1.049891 0.1896514 0.345301 77 : < 1.399105e-05 > ( 21.60612 ) 1.049891 0.1896514 0.345301 78 : < 1.399105e-05 > ( 21.60612 ) 1.049891 0.1896514 0.345301 79 : < 1.313885e-05 > ( 21.60612 ) 1.049891 0.1896514 0.345301 80 : < 1.313885e-05 > ( 21.60612 ) 1.049891 0.1896514 0.345301 81 : < 9.167435e-06 > ( 21.60612 ) 1.049891 0.1896514 0.345301 82 : < 9.167435e-06 > ( 21.60612 ) 1.049891 0.1896514 0.345301 83 : < 9.167435e-06 > ( 21.60612 ) 1.049891 0.1896514 0.345301 84 : < 1.399419e-05 > ( 16.63388 ) 0.9720682 0.2055296 0.3462078 85 : < 1.399419e-05 > ( 16.63388 ) 0.9720682 0.2055296 0.3462078 86 : < 1.399419e-05 > ( 16.63388 ) 0.9720682 0.2055296 0.3462078 87 : < 1.399419e-05 > ( 16.63388 ) 0.9720682 0.2055296 0.3462078 88 : < 1.399419e-05 > ( 16.63388 ) 0.9720682 0.2055296 0.3462078 89 : < 1.399419e-05 > ( 16.63388 ) 0.9720682 0.2055296 0.3462078 90 : < 1.399419e-05 > ( 16.63388 ) 0.9720682 0.2055296 0.3462078 91 : < 1.399419e-05 > ( 16.63388 ) 0.9720682 0.2055296 0.3462078 92 : < 1.399419e-05 > ( 16.63388 ) 0.9720682 0.2055296 0.3462078 93 : < 1.399419e-05 > ( 16.63388 ) 0.9720682 0.2055296 0.3462078 94 : < 1.399419e-05 > ( 16.63388 ) 0.9720682 0.2055296 0.3462078 95 : < 1.399419e-05 > ( 16.63388 ) 0.9720682 0.2055296 0.3462078 96 : < 1.054025e-05 > ( 16.63388 ) 0.9720682 0.2055296 0.3462078 97 : < 9.901587e-06 > ( 16.63388 ) 0.9720682 0.2055296 0.3462078 98 : < 8.687071e-06 > ( 16.63388 ) 0.9720682 0.2055296 0.3462078 99 : < 8.365209e-06 > ( 16.63388 ) 0.9720682 0.2055296 0.3462078 100 : < 8.220149e-06 > ( 16.63388 ) 0.9720682 0.2055296 0.3462078 101 : < 6.218003e-06 > ( 14.77344 ) 0.8998556 0.2261467 0.3049462 102 : < 5.143315e-06 > ( 14.77344 ) 0.8998556 0.2261467 0.3049462 103 : < 5.143315e-06 > ( 14.77344 ) 0.8998556 0.2261467 0.3049462 104 : < 5.143315e-06 > ( 14.77344 ) 0.8998556 0.2261467 0.3049462 105 : < 3.398403e-06 > ( 14.33543 ) 0.8989068 0.227165 0.3331918 106 : < 3.451122e-06 > ( 14.22877 ) 0.8989068 0.227165 0.3389508 107 : < 3.451122e-06 > ( 14.22877 ) 0.8989068 0.227165 0.3389508 108 : < 3.105081e-06 > ( 14.22877 ) 0.8989068 0.227165 0.3389508 109 : < 3.413998e-06 > ( 13.91054 ) 0.8748452 0.2296175 0.3189887 110 : < 2.37363e-06 > ( 13.74816 ) 0.8947597 0.227165 0.3389508 111 : < 2.37363e-06 > ( 13.74816 ) 0.8947597 0.227165 0.3389508 112 : < 2.37363e-06 > ( 13.74816 ) 0.8947597 0.227165 0.3389508 113 : < 2.104666e-06 > ( 13.74816 ) 0.8947597 0.227165 0.3389508 114 : < 1.683264e-06 > ( 13.74816 ) 0.8947597 0.227165 0.3389508 115 : < 1.035148e-06 > ( 13.72196 ) 0.8921428 0.2276308 0.3387805 116 : < 9.497853e-07 > ( 13.69314 ) 0.8910717 0.2275884 0.3347308 117 : < 8.914318e-07 > ( 13.69314 ) 0.8910717 0.2275884 0.3347308 118 : < 8.78976e-07 > ( 13.69314 ) 0.8910717 0.2275884 0.3347308 119 : < 6.091806e-07 > ( 13.69314 ) 0.8910717 0.2275884 0.3347308 120 : < 3.594059e-07 > ( 13.69314 ) 0.8910717 0.2275884 0.3347308 121 : < 3.625165e-07 > ( 13.68993 ) 0.8926421 0.2275628 0.3424516 122 : < 3.625165e-07 > ( 13.68993 ) 0.8926421 0.2275628 0.3424516 123 : < 2.172411e-07 > ( 13.68993 ) 0.8926421 0.2275628 0.3424516 124 : < 2.181822e-07 > ( 13.67633 ) 0.8692439 0.2314282 0.3526782 125 : < 1.701219e-07 > ( 13.64665 ) 0.8830578 0.2289018 0.333503 126 : < 1.108531e-07 > ( 13.64665 ) 0.8830578 0.2289018 0.333503 127 : < 1.594257e-07 > ( 13.59661 ) 0.8874269 0.2282339 0.3396973 128 : < 1.594257e-07 > ( 13.59661 ) 0.8874269 0.2282339 0.3396973 129 : < 1.498057e-07 > ( 13.59661 ) 0.8874269 0.2282339 0.3396973 130 : < 1.498057e-07 > ( 13.59661 ) 0.8874269 0.2282339 0.3396973 131 : < 1.498057e-07 > ( 13.59661 ) 0.8874269 0.2282339 0.3396973 132 : < 1.460906e-07 > ( 13.59661 ) 0.8874269 0.2282339 0.3396973 133 : < 1.460906e-07 > ( 13.59661 ) 0.8874269 0.2282339 0.3396973 134 : < 1.776607e-07 > ( 13.53191 ) 0.8755415 0.2301088 0.3437866 135 : < 1.530891e-07 > ( 13.53191 ) 0.8755415 0.2301088 0.3437866 136 : < 1.489112e-07 > ( 13.53191 ) 0.8755415 0.2301088 0.3437866 137 : < 1.489112e-07 > ( 13.53191 ) 0.8755415 0.2301088 0.3437866 138 : < 9.388662e-08 > ( 13.53191 ) 0.8755415 0.2301088 0.3437866 > showProc.time() Time (user system elapsed): 1.06 0 1.06 > > plot(y ~ x, d.exp.Hlev, main = "Data = d.exp.Hlev") > cE <- curve(Expo(x, a=1, b=0.2), 0, par("usr")[2], n=1+2^9, col=cTr, lwd=2, lty=2, add=TRUE) > x.H <- seq(par("usr")[1], par("usr")[2], length.out = 256) > ll <- length(m1 <- sapply(ls.str(patt="^Rf.Hlev"), get, simplify=FALSE)) > .tmp <- lapply(m1, function(.) lines(x.H, predict(., list(x=x.H)))) > lines(x.H, predict(Cfit.no.Hlev, list(x=x.H)), col=cLS, lwd=3)## L.S.() > if(Cfit.Hlev..ok) { + lines(x.H, predict(Cfit.Hlev, list(x=x.H)), col=cLS, lty=2)## L.S. + legend("topleft", c("true", "LS [w/o outl]", "LS", names(m1)), + lwd=c(2,3, rep(1,1+ll)), lty=c(2,1,2, rep(1,ll)), + col=c(cTr, cLS,cLS, rep(par("fg"),ll)), bty="n", inset=.01) + } else { + cat("no Cfit.Hlev lines as nls() failed there\n") + cat(" : legend(...) !?\n") + } > showProc.time() Time (user system elapsed): 0.02 0 0.02 > > cfcl <- coef(Cfit) > if(Cfit.no..ok) { + cfcl.n.o <- coef(Cfit.no.out) + } else { cfcl.n.o <- cfcl } # use substitute - for code below > cfcl.n.H <- coef(Cfit.no.Hlev) > ## no outliers present > assert.EQ(coef(Rfit.MM.S.bisquare), cfcl, tol = 0.01, giveRE=TRUE) Mean relative difference: 0.003040325 > if(doExtras) { + assert.EQ(coef(Rfit.MM.S.lqq), cfcl, tol = 0.01, giveRE=TRUE) + assert.EQ(coef(Rfit.MM.S.optimal), cfcl, tol = 0.01, giveRE=TRUE) + assert.EQ(coef(Rfit.MM.S.hampel), cfcl, tol = 0.01, giveRE=TRUE) + } > assert.EQ(coef(Rfit.MM.lts.bisquare), cfcl, tol = 0.01, giveRE=TRUE) Mean relative difference: 0.003004411 > assert.EQ(coef(Rfit.MM.lts.lqq), cfcl, tol = 0.01, giveRE=TRUE) Mean relative difference: 0.0004708759 > assert.EQ(coef(Rfit.MM.lts.optimal), cfcl, tol = 0.01, giveRE=TRUE) Mean relative difference: 0.005092942 > assert.EQ(coef(Rfit.MM.lts.hampel), cfcl, tol = 0.01, giveRE=TRUE) Mean relative difference: 0.002536892 > assert.EQ(coef(Rfit.tau.bisquare), cfcl, tol = 0.02, giveRE=TRUE)# 0.009873 Mean relative difference: 0.01018781 > assert.EQ(coef(Rfit.tau.optimal), cfcl, tol = 0.01, giveRE=TRUE) Mean relative difference: 0.00014252 > assert.EQ(coef(Rfit.CM)[-3], cfcl, tol = 0.01, giveRE=TRUE) Mean relative difference: 0.002979679 > assert.EQ(coef(Rfit.mtl)[-3], cfcl, tol = 0.02, giveRE=TRUE) Mean relative difference: 0.002798815 > ## 40% outliers present -- compare with L.S.(good.data) > if(doExtras) { + assert.EQ(coef(Rf.out.MM.S.bisquare), cfcl.n.o, tol = 7e-4, giveRE=TRUE) + assert.EQ(coef(Rf.out.MM.S.lqq), cfcl.n.o, tol = 1e-5, giveRE=TRUE) + assert.EQ(coef(Rf.out.MM.S.optimal), cfcl.n.o, tol = 1e-5, giveRE=TRUE) + assert.EQ(coef(Rf.out.MM.S.hampel), cfcl.n.o, tol = 1e-5, giveRE=TRUE) + } > assert.EQ(coef(Rf.out.MM.lts.bisquare), cfcl.n.o, tol = 6e-4, giveRE=TRUE) Mean relative difference: 0.0005948213 > assert.EQ(coef(Rf.out.MM.lts.lqq), cfcl.n.o, tol = 1e-5, giveRE=TRUE) Mean relative difference: 7.09538e-06 > assert.EQ(coef(Rf.out.MM.lts.optimal), cfcl.n.o, tol = 1e-5, giveRE=TRUE) Mean relative difference: 7.658257e-06 > assert.EQ(coef(Rf.out.MM.lts.hampel), cfcl.n.o, tol = 1e-5, giveRE=TRUE) Mean relative difference: 7.683842e-06 > assert.EQ(coef(Rf.out.tau.bisquare), cfcl.n.o, tol = .007, giveRE=TRUE) Mean relative difference: 0.005147822 > assert.EQ(coef(Rf.out.tau.optimal), cfcl.n.o, tol = .002, giveRE=TRUE) Mean relative difference: 0.001682991 > assert.EQ(coef(Rf.out.CM)[-3], cfcl.n.o, tol = .012, giveRE=TRUE)# 0.00708,0.01079 Mean relative difference: 0.00881462 > assert.EQ(coef(Rf.out.mtl)[-3], cfcl.n.o, tol = .002, giveRE=TRUE)# better in 64b Mean relative difference: 0.001445415 > ## presence of high leverage point outliers -- compare with LS(good.data) > if(doExtras) { + assert.EQ(coef(Rf.Hlev.MM.S.bisquare), cfcl.n.H, tol = .01, giveRE=TRUE) + assert.EQ(coef(Rf.Hlev.MM.S.lqq), cfcl.n.H, tol = .02, giveRE=TRUE) + assert.EQ(coef(Rf.Hlev.MM.S.optimal), cfcl.n.H, tol = .005, giveRE=TRUE) + assert.EQ(coef(Rf.Hlev.MM.S.hampel), cfcl.n.H, tol = .02, giveRE=TRUE) + } > assert.EQ(coef(Rf.Hlev.MM.lts.bisquare),cfcl.n.H, tol = .01, giveRE=TRUE) Mean relative difference: 0.007644998 > assert.EQ(coef(Rf.Hlev.MM.lts.lqq), cfcl.n.H, tol = .015, giveRE=TRUE) Mean relative difference: 0.01077487 > assert.EQ(coef(Rf.Hlev.MM.lts.optimal), cfcl.n.H, tol = .002, giveRE=TRUE) Mean relative difference: 0.001024011 > assert.EQ(coef(Rf.Hlev.MM.lts.hampel), cfcl.n.H, tol = .02, giveRE=TRUE) Mean relative difference: 0.01395353 > assert.EQ(coef(Rf.Hlev.tau.bisquare), cfcl.n.H, tol = .05, giveRE=TRUE)# 0.0363, 0.0415 Mean relative difference: 0.02634211 > assert.EQ(coef(Rf.Hlev.tau.optimal), cfcl.n.H, tol = .03, giveRE=TRUE) Mean relative difference: 0.01064765 > assert.EQ(coef(Rf.Hlev.CM)[-3], cfcl.n.H, tol = .12, giveRE=TRUE)# 0.032, 0.082 Mean relative difference: 0.009277506 > assert.EQ(coef(Rf.Hlev.mtl)[-3], cfcl.n.H, tol = .08, giveRE=TRUE) Mean relative difference: 0.05156692 > > length(mods <- sapply(ls.str(patt="^Rf"), get, simplify=FALSE)) # 36 [1] 25 > is.conv <- sapply(mods, `[[`, "status") == "converged" > prblm <- mods[!is.conv] > if(length(prblm)) { + cat("\n*** NON-converged model fits:\n") + print(prblm) + mods <- mods[is.conv] + } else cat("\n All models converged\n") All models converged > > ## Now, all mods are converged ----------- > > dKnd <- as.factor(vapply(mods, function(.m.) + as.character(getCall(.m.)[["data"]]), "")) > table(dKnd) ## dKnd d.exp.Hlev d.exp30 d.exp40out 8 9 8 > (iKnd <- setNames(seq_len(nlevels(dKnd)), levels(dKnd))) d.exp.Hlev d.exp30 d.exp40out 1 2 3 > > ## Coefficients: Some have 'sigma', some not: > pcf <- vapply(lcf <- lapply(mods, coef), length, 1) > table(pcf) ## 2 and 3 pcf 2 3 19 6 > stopifnot(min(pcf) + 1 == max(pcf)) # +1 : those which have 'sigma > pp <- min(pcf) > ccf <- t(simplify2array(lapply(lcf, `[`, 1:max(pcf)))) > ## take the "Scale" for those that do not have 'sigma' among coef(): > i.n <- is.na(ccf[,"sigma"]) > ccf[i.n, "sigma"] <- vapply(mods[i.n], `[[`, 0, "Scale") > ## not yet: vapply(mods[i.n], sigma, 0.) > ccf a b sigma Rf.Hlev.CM 0.9120862 0.2201427 3.2251456 Rf.Hlev.MM.lts.bisquare 0.9119031 0.2180941 0.9104102 Rf.Hlev.MM.lts.hampel 0.9056888 0.2189335 0.9226696 Rf.Hlev.MM.lts.lqq 0.9087898 0.2184883 0.9265427 Rf.Hlev.MM.lts.optimal 0.9182742 0.2169890 0.9059446 Rf.Hlev.mtl 0.8755415 0.2301088 0.3437866 Rf.Hlev.tau.bisquare 0.9476494 0.2145860 1.2734884 Rf.Hlev.tau.optimal 0.9311202 0.2164535 1.0953166 Rf.out.CM 0.9741043 0.2001078 3.6869833 Rf.out.MM.lts.bisquare 0.9828761 0.1992326 2.4043186 Rf.out.MM.lts.hampel 0.9835150 0.1991592 2.5777817 Rf.out.MM.lts.lqq 0.9835144 0.1991593 2.5428862 Rf.out.MM.lts.optimal 0.9835150 0.1991592 2.2454196 Rf.out.mtl 0.9819907 0.1993515 0.7891589 Rf.out.tau.bisquare 0.9780900 0.1998069 5.7465801 Rf.out.tau.optimal 0.9817293 0.1993704 4.1313023 Rfit.CM 0.9995511 0.2005112 3.8722961 Rfit.MM.S.bisquare 0.9995623 0.2004497 0.7503993 Rfit.MM.lts.bisquare 0.9995257 0.2004563 0.7535988 Rfit.MM.lts.hampel 0.9991106 0.2006029 0.7635464 Rfit.MM.lts.lqq 0.9969093 0.2008812 0.7613701 Rfit.MM.lts.optimal 0.9911139 0.2017249 0.7420144 Rfit.mtl 0.9987386 0.1999197 0.2792836 Rfit.tau.bisquare 1.0069796 0.1992269 0.8105042 Rfit.tau.optimal 0.9963417 0.2010482 0.8332127 > ## well, the 'sigma's are definitely *not* unbiased estimates of > ## true sqrt(var(eps)) ... [FIXME] > ## --> indeed, this can be found in the CM paper [TODO: write more here] > > plot(ccf[,1:2], pch = as.integer(dKnd))## use 'method' to get color > legend("topright", inset=.01, names(iKnd), pch = iKnd) > points(rbind(cfcl.n.H, cfcl, cfcl.n.o), # <- order from iKind + col=adjustcolor("tomato",.5), cex=2, pch=1:3, lwd=5) > ## optional > labs <- sub("^Rfit\\.", '', sub("^Rf\\.[A-Za-z]+\\.", '', rownames(ccf))) > labs <- sub("hampel$", "Ham", sub("optimal$", "opt", sub("bisquare$", "biS", labs))) > labs [1] "CM" "MM.lts.biS" "MM.lts.Ham" "MM.lts.lqq" "MM.lts.opt" [6] "mtl" "tau.biS" "tau.opt" "CM" "MM.lts.biS" [11] "MM.lts.Ham" "MM.lts.lqq" "MM.lts.opt" "mtl" "tau.biS" [16] "tau.opt" "CM" "MM.S.biS" "MM.lts.biS" "MM.lts.Ham" [21] "MM.lts.lqq" "MM.lts.opt" "mtl" "tau.biS" "tau.opt" > text(ccf[,1:2], labs, cex=0.75, col=adjustcolor(1, 0.5), + adj= -1/5, srt=75, xpd=NA) > points(rbind(cfcl), col=adjustcolor("tomato",.5), cex=2, pch=3, lwd=5) > showProc.time() Time (user system elapsed): 0.02 0 0.01 > > > ###------- Extended Tests for the DNase1 example from >>>> ../man/nlrob-algos.Rd <<<< > ### ===================== > DNase1 <- DNase[DNase$Run == 1,] > form <- density ~ Asym/(1 + exp(( xmid -log(conc) )/scal )) > pnms <- c("Asym", "xmid", "scal") > psNms <- c(pnms, "sigma") > > ##' a version that recycles x: > setNames. <- function(x, nm) setNames(rep_len(x, length(nm)), nm) > > ## for comparisons, later: > all.eq.mod <- function(m1, m2, sub=FALSE, excl = c("call", "ctrl"), ...) { + nm1 <- names(m1) + stopifnot(if(sub) nm1 %in% names(m2) else nm1 == names(m2)) + ni <- if(sub) + nm1[is.na(match(nm1, c("call","ctrl")))] + else is.na(match(names(m1), excl))## <<- all but those with names in 'excl' + all.equal(m1[ni], m2[ni], ...) + } > > set.seed(47) # as these by default use randomized optimization: > fMM <- robustbase:::nlrob.MM(form, data = DNase1, + lower = setNames.(0, pnms), upper = 3, + ## call to nlrob.control to pass 'optim.control': + ctrl = nlrob.control("MM", optim.control = list(trace = 1), + optArgs = list(trace = TRUE))) 1 : < 0.06568043 > ( 0.030657 ) 2.699558 1.662672 1.071085 2 : < 0.05135563 > ( 0.01986514 ) 2.367089 1.520246 1.071085 3 : < 0.03176846 > ( 0.01986514 ) 2.367089 1.520246 1.071085 4 : < 0.03073827 > ( 0.01986514 ) 2.367089 1.520246 1.071085 5 : < 0.02446155 > ( 0.01986514 ) 2.367089 1.520246 1.071085 6 : < 0.02158731 > ( 0.01986514 ) 2.367089 1.520246 1.071085 7 : < 0.02076038 > ( 0.01986514 ) 2.367089 1.520246 1.071085 8 : < 0.01409694 > ( 0.01986514 ) 2.367089 1.520246 1.071085 9 : < 0.01014985 > ( 0.01985212 ) 2.303141 1.388635 1.028954 10 : < 0.009128672 > ( 0.01985212 ) 2.303141 1.388635 1.028954 11 : < 0.008785776 > ( 0.01985212 ) 2.303141 1.388635 1.028954 12 : < 0.008176153 > ( 0.01882007 ) 1.429895 0.7101177 0.9320048 13 : < 0.007138518 > ( 0.01882007 ) 1.429895 0.7101177 0.9320048 14 : < 0.00663419 > ( 0.01694984 ) 2.264548 1.373464 1.007218 15 : < 0.004990433 > ( 0.01218113 ) 2.634336 1.786663 1.130746 16 : < 0.004450077 > ( 0.01218113 ) 2.634336 1.786663 1.130746 17 : < 0.003781576 > ( 0.01218113 ) 2.634336 1.786663 1.130746 18 : < 0.003439636 > ( 0.01218113 ) 2.634336 1.786663 1.130746 19 : < 0.002096688 > ( 0.01218113 ) 2.634336 1.786663 1.130746 20 : < 0.001858781 > ( 0.01218113 ) 2.634336 1.786663 1.130746 21 : < 0.001466238 > ( 0.01218113 ) 2.634336 1.786663 1.130746 22 : < 0.001418615 > ( 0.01218113 ) 2.634336 1.786663 1.130746 23 : < 0.001301314 > ( 0.01218113 ) 2.634336 1.786663 1.130746 24 : < 0.001301314 > ( 0.01218113 ) 2.634336 1.786663 1.130746 25 : < 0.001223393 > ( 0.01218113 ) 2.634336 1.786663 1.130746 26 : < 0.0009908125 > ( 0.01218113 ) 2.634336 1.786663 1.130746 27 : < 0.000939791 > ( 0.01218113 ) 2.634336 1.786663 1.130746 28 : < 0.0009356657 > ( 0.01218113 ) 2.634336 1.786663 1.130746 29 : < 0.0009316897 > ( 0.01218113 ) 2.634336 1.786663 1.130746 30 : < 0.0008635449 > ( 0.01211585 ) 2.594124 1.780671 1.127007 31 : < 0.0009441278 > ( 0.01153015 ) 2.574215 1.73965 1.117878 32 : < 0.0009441278 > ( 0.01153015 ) 2.574215 1.73965 1.117878 33 : < 0.0009234668 > ( 0.01153015 ) 2.574215 1.73965 1.117878 34 : < 0.001017899 > ( 0.0104385 ) 2.555245 1.744374 1.124387 35 : < 0.0009771945 > ( 0.01011045 ) 2.715392 1.880491 1.15727 36 : < 0.001117598 > ( 0.00935059 ) 2.720104 1.882524 1.155785 37 : < 0.0008365975 > ( 0.00935059 ) 2.720104 1.882524 1.155785 38 : < 0.0007486373 > ( 0.00935059 ) 2.720104 1.882524 1.155785 39 : < 0.0007486373 > ( 0.00935059 ) 2.720104 1.882524 1.155785 40 : < 0.0006539709 > ( 0.00935059 ) 2.720104 1.882524 1.155785 41 : < 0.0006052928 > ( 0.00935059 ) 2.720104 1.882524 1.155785 42 : < 0.0008283781 > ( 0.007374737 ) 2.70406 1.863656 1.151602 43 : < 0.0007742664 > ( 0.007374737 ) 2.70406 1.863656 1.151602 44 : < 0.0006118239 > ( 0.007374737 ) 2.70406 1.863656 1.151602 45 : < 0.000559146 > ( 0.007374737 ) 2.70406 1.863656 1.151602 46 : < 0.0005460518 > ( 0.00735965 ) 2.69992 1.859164 1.149771 47 : < 0.0005050816 > ( 0.00735965 ) 2.69992 1.859164 1.149771 48 : < 0.0004511749 > ( 0.00735965 ) 2.69992 1.859164 1.149771 49 : < 0.0004511749 > ( 0.00735965 ) 2.69992 1.859164 1.149771 50 : < 0.0004308587 > ( 0.00735965 ) 2.69992 1.859164 1.149771 51 : < 0.0004021508 > ( 0.007323208 ) 2.705071 1.86404 1.151824 52 : < 0.0004021508 > ( 0.007323208 ) 2.705071 1.86404 1.151824 53 : < 0.00032253 > ( 0.007323208 ) 2.705071 1.86404 1.151824 54 : < 0.0002989224 > ( 0.007268853 ) 2.701306 1.860916 1.150708 55 : < 0.0002612988 > ( 0.007268853 ) 2.701306 1.860916 1.150708 56 : < 0.0002094852 > ( 0.007268853 ) 2.701306 1.860916 1.150708 57 : < 0.0002141296 > ( 0.007162363 ) 2.688279 1.849162 1.147415 58 : < 0.0001640686 > ( 0.007162363 ) 2.688279 1.849162 1.147415 59 : < 0.0001162773 > ( 0.007162363 ) 2.688279 1.849162 1.147415 60 : < 0.0001162773 > ( 0.007162363 ) 2.688279 1.849162 1.147415 61 : < 0.0001162773 > ( 0.007162363 ) 2.688279 1.849162 1.147415 62 : < 7.936557e-05 > ( 0.007138531 ) 2.695107 1.85494 1.149358 63 : < 4.461414e-05 > ( 0.007138531 ) 2.695107 1.85494 1.149358 64 : < 4.289425e-05 > ( 0.007138531 ) 2.695107 1.85494 1.149358 65 : < 2.867267e-05 > ( 0.007138531 ) 2.695107 1.85494 1.149358 66 : < 2.209151e-05 > ( 0.007138531 ) 2.695107 1.85494 1.149358 67 : < 1.419382e-05 > ( 0.007138531 ) 2.695107 1.85494 1.149358 68 : < 1.617992e-05 > ( 0.007115951 ) 2.693241 1.853284 1.148951 69 : < 1.431852e-05 > ( 0.007115951 ) 2.693241 1.853284 1.148951 70 : < 1.173847e-05 > ( 0.007114441 ) 2.687203 1.847581 1.147502 71 : < 1.153084e-05 > ( 0.007114441 ) 2.687203 1.847581 1.147502 72 : < 1.112424e-05 > ( 0.007114441 ) 2.687203 1.847581 1.147502 73 : < 8.432844e-06 > ( 0.007114441 ) 2.687203 1.847581 1.147502 74 : < 6.626809e-06 > ( 0.007108962 ) 2.690307 1.850949 1.148608 75 : < 5.855829e-06 > ( 0.007108962 ) 2.690307 1.850949 1.148608 76 : < 6.557656e-06 > ( 0.007099157 ) 2.692359 1.852365 1.148664 77 : < 6.341923e-06 > ( 0.007099157 ) 2.692359 1.852365 1.148664 78 : < 4.799942e-06 > ( 0.007099157 ) 2.692359 1.852365 1.148664 79 : < 3.70222e-06 > ( 0.007099157 ) 2.692359 1.852365 1.148664 80 : < 3.299484e-06 > ( 0.007099157 ) 2.692359 1.852365 1.148664 81 : < 2.867961e-06 > ( 0.007095899 ) 2.691273 1.851573 1.148496 82 : < 2.419718e-06 > ( 0.007095899 ) 2.691273 1.851573 1.148496 83 : < 3.092519e-06 > ( 0.007091715 ) 2.691252 1.851506 1.148814 84 : < 3.092519e-06 > ( 0.007091715 ) 2.691252 1.851506 1.148814 85 : < 2.59233e-06 > ( 0.007091715 ) 2.691252 1.851506 1.148814 86 : < 2.494284e-06 > ( 0.007091715 ) 2.691252 1.851506 1.148814 87 : < 2.494284e-06 > ( 0.007091715 ) 2.691252 1.851506 1.148814 88 : < 2.494284e-06 > ( 0.007091715 ) 2.691252 1.851506 1.148814 89 : < 1.951421e-06 > ( 0.007090575 ) 2.691043 1.851372 1.148696 90 : < 1.951421e-06 > ( 0.007090575 ) 2.691043 1.851372 1.148696 91 : < 1.469318e-06 > ( 0.007090575 ) 2.691043 1.851372 1.148696 92 : < 1.038646e-06 > ( 0.007090575 ) 2.691043 1.851372 1.148696 93 : < 8.753423e-07 > ( 0.007090491 ) 2.691199 1.851383 1.148837 iter 10 value 732.941462 final value 732.397073 converged > showProc.time() Time (user system elapsed): 0.9 0.01 0.93 > > if(doExtras) { + ftau <- robustbase:::nlrob.tau(form, data = DNase1, + lower= setNames.(0, pnms), upper= 3, trace=TRUE) + ## + fCM <- robustbase:::nlrob.CM (form, data = DNase1, + lower= setNames.(0, psNms), upper= 3, trace=TRUE) + ## + fmtl <- robustbase:::nlrob.mtl(form, data = DNase1, + lower= setNames.(0, psNms), upper= 3, trace=TRUE) + ## + mods <- list(MM=fMM, tau=ftau, CM=fCM, MTL=fmtl) + print(sts <- sapply(mods, `[[`, "status")) + stopifnot(sts == "converged") + + print(sapply(mods, `[[`, "data")) # currently 'language' %% FIXME + + print(sapply(mods, `[[`, "coefficients")) # nice matrix + showProc.time() + } > ## Compare with traditional M-estimate, a) started robustly b) psi = Tukey's: > fM <- nlrob(formula(fMM), data=eval(fMM$data), start = coef(fMM), + psi = .Mwgt.psi1("bisquare"), trace = TRUE) robust iteration 1 0.001151704 (1.86e-01): par = (2.561138 1.730322 1.112787) 0.001113843 (9.35e-03): par = (2.528926 1.694972 1.101376) 0.001113745 (6.63e-05): par = (2.530272 1.696145 1.101659) 0.001113745 (1.07e-06): par = (2.53026 1.696132 1.101654) --> irls.delta(previous, resid) = 0.0890318 -- *not* converged robust iteration 2 0.001403458 (1.41e-01): par = (2.53026 1.696132 1.101654) 0.001376449 (5.85e-03): par = (2.505189 1.668215 1.093201) 0.001376402 (3.91e-05): par = (2.506056 1.668969 1.093384) 0.001376402 (6.00e-07): par = (2.506048 1.668961 1.093381) --> irls.delta(previous, resid) = 0.0782933 -- *not* converged robust iteration 3 0.001771648 (1.47e-01): par = (2.506048 1.668961 1.093381) 0.001734285 (6.92e-03): par = (2.47793 1.637687 1.084693) 0.001734202 (2.19e-05): par = (2.478792 1.638387 1.084841) 0.001734202 (2.52e-07): par = (2.478787 1.638382 1.084839) --> irls.delta(previous, resid) = 0.0942244 -- *not* converged robust iteration 4 0.002327706 (1.70e-01): par = (2.478787 1.638382 1.084839) 0.002262441 (9.12e-03): par = (2.444947 1.600353 1.074457) 0.002262253 (8.94e-06): par = (2.445935 1.601097 1.074576) --> irls.delta(previous, resid) = 0.125083 -- *not* converged robust iteration 5 0.002375777 (7.00e-02): par = (2.445935 1.601097 1.074576) 0.002364223 (1.35e-03): par = (2.43303 1.58614 1.070358) 0.002364219 (9.27e-07): par = (2.433202 1.586277 1.070384) --> irls.delta(previous, resid) = 0.0556974 -- *not* converged robust iteration 6 0.002562082 (6.45e-02): par = (2.433202 1.586277 1.070384) 0.002551476 (1.17e-03): par = (2.421044 1.57219 1.066458) 0.002551472 (1.99e-07): par = (2.421173 1.572288 1.066473) --> irls.delta(previous, resid) = 0.0537787 -- *not* converged robust iteration 7 0.002952554 (9.54e-02): par = (2.421173 1.572288 1.066473) 0.002925917 (2.57e-03): par = (2.402659 1.550824 1.060525) 0.002925898 (1.84e-06): par = (2.402888 1.550976 1.060531) --> irls.delta(previous, resid) = 0.0838695 -- *not* converged robust iteration 8 0.003345785 (9.84e-02): par = (2.402888 1.550976 1.060531) 0.003313638 (2.58e-03): par = (2.383987 1.528804 1.054326) 0.003313616 (1.72e-06): par = (2.384174 1.528911 1.054311) --> irls.delta(previous, resid) = 0.0910549 -- *not* converged robust iteration 9 0.003557050 (6.62e-02): par = (2.384174 1.528911 1.054311) 0.003541462 (1.09e-03): par = (2.371903 1.514284 1.050158) 0.003541457 (6.01e-07): par = (2.37194 1.514285 1.050134) --> irls.delta(previous, resid) = 0.0632662 -- *not* converged robust iteration 10 0.003657973 (3.61e-02): par = (2.37194 1.514285 1.050134) 0.003653196 (3.34e-04): par = (2.365452 1.506462 1.047891) 0.003653196 (1.05e-06): par = (2.365434 1.506433 1.047873) --> irls.delta(previous, resid) = 0.0349519 -- *not* converged robust iteration 11 0.003681549 (1.58e-02): par = (2.365434 1.506433 1.047873) 0.003680628 (9.50e-05): par = (2.362727 1.503117 1.046909) 0.003680628 (6.44e-07): par = (2.36271 1.503097 1.0469) --> irls.delta(previous, resid) = 0.0153554 -- *not* converged robust iteration 12 0.003691710 (6.80e-03): par = (2.36271 1.503097 1.0469) 0.003691539 (3.85e-05): par = (2.361561 1.501683 1.046487) 0.003691539 (3.38e-07): par = (2.361552 1.501674 1.046483) --> irls.delta(previous, resid) = 0.00662493 -- *not* converged robust iteration 13 0.003696077 (2.91e-03): par = (2.361552 1.501674 1.046483) 0.003696046 (1.66e-05): par = (2.361064 1.501071 1.046307) 0.003696046 (1.04e-07): par = (2.36106 1.501066 1.046305) --> irls.delta(previous, resid) = 0.00283586 -- *not* converged robust iteration 14 0.003697945 (1.24e-03): par = (2.36106 1.501066 1.046305) 0.003697939 (7.19e-06): par = (2.360852 1.50081 1.04623) --> irls.delta(previous, resid) = 0.00120506 -- *not* converged robust iteration 15 0.003698734 (5.31e-04): par = (2.360852 1.50081 1.04623) 0.003698733 (3.10e-06): par = (2.360762 1.500699 1.046197) --> irls.delta(previous, resid) = 0.000515829 -- *not* converged robust iteration 16 0.003699076 (2.27e-04): par = (2.360762 1.500699 1.046197) 0.003699076 (1.42e-06): par = (2.360723 1.500652 1.046183) --> irls.delta(previous, resid) = 0.000220909 -- *not* converged robust iteration 17 0.003699222 (9.74e-05): par = (2.360723 1.500652 1.046183) 0.003699222 (4.75e-07): par = (2.360707 1.500631 1.046177) --> irls.delta(previous, resid) = 9.46532e-05 -- *not* converged robust iteration 18 0.003699285 (4.17e-05): par = (2.360707 1.500631 1.046177) 0.003699285 (2.57e-07): par = (2.3607 1.500623 1.046175) --> irls.delta(previous, resid) = 4.04661e-05 -- *not* converged robust iteration 19 0.003699312 (1.78e-05): par = (2.3607 1.500623 1.046175) 0.003699312 (1.83e-07): par = (2.360697 1.500619 1.046174) --> irls.delta(previous, resid) = 1.73262e-05 -- *not* converged robust iteration 20 0.003699324 (7.69e-06): par = (2.360697 1.500619 1.046174) > rbind(M=coef(fM), MM=coef(fMM)) # "similar" ... well, no: the sigma's get much different Asym xmid scal M 2.360697 1.500619 1.046174 MM 2.561138 1.730322 1.112787 > ## stopifnot(%%____FIXME___ > all.equal(coef(fM), coef(fMM), tolerance = 1e-4) [1] "Mean relative difference: 0.1012243" > ## ) # had 3.26e-5 > ## FIXME: nlrob( "M") should allow to keep specify an initial sigma *and* keep that fixed > showProc.time() Time (user system elapsed): 0.06 0 0.06 > > if(doExtras) { + ### Now call the above methods via nlrob(): + set.seed(47) # (same as above) + ## without "sigma" + gMM <- nlrob(form, data = DNase1, method = "MM", + lower = setNames(c(0,0,0), pnms), upper = 3) + gtau <- nlrob(form, data = DNase1, method = "tau", + lower = setNames(c(0,0,0), pnms), upper = 3) + ## those with "sigma" -> must be in (lower, upper), too : + gCM <- nlrob(form, data = DNase1, method = "CM", + lower = setNames(c(0,0,0,0), psNms), upper = 3) + gmtl <- nlrob(form, data = DNase1, method = "mtl", + lower = setNames(c(0,0,0,0), psNms), upper = 3) + showProc.time() + + ## list {and test print() for these}: + (mod2 <- list(MM=gMM, tau=gtau, CM=gCM, MTL=gmtl)) + + stopifnot(mapply(all.eq.mod, mods, mod2, sub=TRUE)) + }# only if(doExtras) > > proc.time() user system elapsed 7.96 0.10 8.06