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. > library(robustbase) > > ## instead of relying on system.file("test-tools-1.R", package="Matrix"): > source(system.file("xtraR/test-tools.R", package = "robustbase")) > > ###>> 1 ------------------- family = poisson ------------------------------------ > > ### very simple model [with outliers] > set.seed(113) > y <- rpois(17, lambda = 4) ## -> target: beta_0 = log(E[Y]) = log(4) = 1.386294 > > y[1:2] <- 99:100 # outliers > y [1] 99 100 7 3 5 3 2 5 3 5 1 3 2 3 4 6 3 > rm1 <- glmrob(y ~ 1, family = poisson, trace = TRUE, + acc = 1e-13) # default is just 1e-4 Initial theta: (Intercept) 2.704121 it | d{theta} | rel.change 1 | -0.32 | 0.117708 2 | -0.37 | 0.154689 3 | -0.36 | 0.177974 4 | -0.2 | 0.122032 5 | -0.035 | 0.0243059 6 | -0.0028 | 0.00199645 7 | -0.00019 | 0.000132804 8 | -1.2e-5 | 8.7184e-06 9 | -8.1e-7 | 5.71848e-07 10 | -5.3e-8 | 3.75058e-08 11 | -3.5e-9 | 2.45989e-09 12 | -2.3e-10 | 1.61336e-10 13 | -1.5e-11 | 1.05815e-11 14 | -9.8e-13 | 6.94032e-13 15 | -6.5e-14 | 4.55518e-14 > ## and check the robustness weights: > assert.EQ(c(0.0287933850640724, 0.0284930623638766, + 0.950239140568007, 0.874115394740014), + local({w <- rm1$w.r; w[ w != 1 ] }), tol = 1e-14) > assert.EQ(coef(rm1), c("(Intercept)" = 1.41710946076738),tol = 1e-14) > > cm1 <- glm (y ~ 1, family = poisson, trace = TRUE) Deviance = 1163.408 Iterations - 1 Deviance = 683.2127 Iterations - 2 Deviance = 613.1911 Iterations - 3 Deviance = 610.3937 Iterations - 4 Deviance = 610.3869 Iterations - 5 Deviance = 610.3869 Iterations - 6 > > rmMT <- glmrob(y ~ 1, family = poisson, trace = TRUE, method="MT") Computing initial estimate with 500 sub samples: .................................................. 50 .................................................. 100 .................................................. 150 .................................................. 200 .................................................. 250 .................................................. 300 .................................................. 350 .................................................. 400 .................................................. 450 .................................................. 500 Optim()izing sumaConPesos() final value 3.622702 converged > (sMT <- summary(rmMT)) Call: glmrob(formula = y ~ 1, family = poisson, method = "MT", trace.lev = TRUE) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 1.2937 0.1205 10.73 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Robustness weights w.r * w.x: 2 observations c(1,2) are outliers with |weight| = 0 ( < 0.0059); The remaining 15 ones are summarized as Min. 1st Qu. Median Mean 3rd Qu. Max. 0.6868 0.9084 0.9399 0.9242 0.9915 0.9922 Number of observations: 17 Fitted by method 'MT' (in 14 iterations) (Dispersion parameter for poisson family taken to be 1) No deviance values available Algorithmic parameters: cw acc 2.1e+00 1.0e-06 nsubm maxit 500 200 > > if(FALSE) # for manual digging: + debug(robustbase:::glmrobMqle) > > allresid <- function(obj, types = c("deviance", "pearson", "working", "response")) + { + sapply(types, residuals, object = obj) + } > > okFit <- function(obj, check.attr=FALSE, ...) { + all.equal(obj$y, obj$fitted.values + residuals(obj, "response"), + check.attributes=check.attr, ...) + } > > ## check validity of several methods simultaneously: > y. <- model.response(model.frame(rm1)) > stopifnot(okFit(cm1), okFit(rm1), y. == y) > > alr.c <- allresid(cm1) > alr.r <- allresid(rm1) > > ## MM --- just for now -- > plot(resid(cm1), resid(rm1), asp=1); abline(0,1, col=2) > plot(resid(cm1,type="pearson"), resid(rm1, type="pearson"), asp=1); abline(0,1, col=2) > plot(resid(cm1,type="working"), resid(rm1, type="working"), asp=1); abline(0,1, col=2) > > ## leave away the outliers -- > cm0 <- glm (y ~ 1, family = poisson, trace = TRUE, subset = -(1:2)) Deviance = 10.73791 Iterations - 1 Deviance = 10.29225 Iterations - 2 Deviance = 10.29144 Iterations - 3 Deviance = 10.29144 Iterations - 4 > plot(resid(cm0), resid(rm1)[-(1:2)], asp=1); abline(0,1, col=2) > plot(resid(cm0,type="pearson"), resid(rm1, type="pearson")[-(1:2)], asp=1); abline(0,1, col=2) > plot(resid(cm0,type="working"), resid(rm1, type="working")[-(1:2)], asp=1); abline(0,1, col=2) > plot(resid(cm0,type="response"), resid(rm1, type="response")[-(1:2)], asp=1); abline(0,1, col=2) > > > ## Use weights (slow convergence !) > w2 <- c(rep(1,8), rep(10,9)) > rm2 <- glmrob(y ~ 1, family = poisson, trace = TRUE, + weights = w2, maxit = 500, acc = 1e-10) # default is just 1e-4 Initial theta: (Intercept) 1.676524 it | d{theta} | rel.change 1 | 0.65 | 0.387101 2 | -0.23 | 0.0992927 3 | -0.2 | 0.0953427 4 | -0.032 | 0.0169877 5 | 0.03 | 0.0159738 6 | -0.028 | 0.0146181 7 | 0.026 | 0.0137313 8 | -0.024 | 0.0126081 9 | 0.022 | 0.0118372 10 | -0.021 | 0.010904 11 | 0.019 | 0.0102371 12 | -0.018 | 0.00945966 13 | 0.017 | 0.00888519 14 | -0.016 | 0.00823613 15 | 0.014 | 0.00774296 16 | -0.014 | 0.00720008 17 | 0.013 | 0.00677795 18 | -0.012 | 0.00632314 19 | 0.011 | 0.00596273 20 | -0.011 | 0.00558115 21 | 0.0099 | 0.00527409 22 | -0.0093 | 0.00495355 23 | 0.0088 | 0.00469244 24 | -0.0083 | 0.00442287 25 | 0.0079 | 0.0042012 26 | -0.0075 | 0.00397426 27 | 0.0071 | 0.00378634 28 | -0.0068 | 0.00359512 29 | 0.0064 | 0.00343603 30 | -0.0062 | 0.00327474 31 | 0.0059 | 0.00314023 32 | -0.0056 | 0.00300408 33 | 0.0054 | 0.00289048 34 | -0.0052 | 0.00277544 35 | 0.005 | 0.00267962 36 | -0.0048 | 0.00257497 37 | 0.0047 | 0.00248707 38 | -0.0045 | 0.0023904 39 | 0.0043 | 0.00230839 40 | -0.0042 | 0.00221905 41 | 0.004 | 0.00214257 42 | -0.0039 | 0.00205998 43 | 0.0037 | 0.00198868 44 | -0.0036 | 0.00191231 45 | 0.0035 | 0.00184586 46 | -0.0033 | 0.00177523 47 | 0.0032 | 0.00171331 48 | -0.0031 | 0.00164796 49 | 0.003 | 0.00159029 50 | -0.0029 | 0.00152982 51 | 0.0028 | 0.00147611 52 | -0.0027 | 0.00142014 53 | 0.0026 | 0.00137014 54 | -0.0025 | 0.00131833 55 | 0.0024 | 0.00127179 56 | -0.0023 | 0.00122381 57 | 0.0022 | 0.0011805 58 | -0.0021 | 0.00113606 59 | 0.0021 | 0.00109576 60 | -0.002 | 0.00105461 61 | 0.0019 | 0.00101712 62 | -0.0018 | 0.00097899 63 | 0.0018 | 0.00094412 64 | -0.0017 | 0.000908793 65 | 0.0016 | 0.000876364 66 | -0.0016 | 0.000843628 67 | 0.0015 | 0.000813473 68 | -0.0015 | 0.000783134 69 | 0.0014 | 0.000755097 70 | -0.0014 | 0.000726977 71 | 0.0013 | 0.000700912 72 | -0.0013 | 0.000674846 73 | 0.0012 | 0.000650617 74 | -0.0012 | 0.000626452 75 | 0.0011 | 0.000603933 76 | -0.0011 | 0.000581528 77 | 0.0011 | 0.000560599 78 | -0.001 | 0.000539824 79 | 0.00098 | 0.000520375 80 | -0.00094 | 0.000501111 81 | 0.00091 | 0.000483039 82 | -0.00087 | 0.000465173 83 | 0.00084 | 0.000448382 84 | -0.00081 | 0.000431813 85 | 0.00078 | 0.000416212 86 | -0.00075 | 0.000400844 87 | 0.00073 | 0.00038635 88 | -0.0007 | 0.000372096 89 | 0.00067 | 0.000358632 90 | -0.00065 | 0.00034541 91 | 0.00063 | 0.000332903 92 | -0.0006 | 0.000320637 93 | 0.00058 | 0.000309019 94 | -0.00056 | 0.000297641 95 | 0.00054 | 0.00028685 96 | -0.00052 | 0.000276293 97 | 0.0005 | 0.000266271 98 | -0.00048 | 0.000256477 99 | 0.00046 | 0.000247169 100 | -0.00045 | 0.000238082 101 | 0.00043 | 0.000229437 102 | -0.00042 | 0.000221006 103 | 0.0004 | 0.000212977 104 | -0.00039 | 0.000205154 105 | 0.00037 | 0.000197699 106 | -0.00036 | 0.00019044 107 | 0.00034 | 0.000183516 108 | -0.00033 | 0.000176781 109 | 0.00032 | 0.000170351 110 | -0.00031 | 0.000164101 111 | 0.0003 | 0.000158131 112 | -0.00029 | 0.000152331 113 | 0.00028 | 0.000146787 114 | -0.00027 | 0.000141405 115 | 0.00026 | 0.000136257 116 | -0.00025 | 0.000131262 117 | 0.00024 | 0.000126483 118 | -0.00023 | 0.000121847 119 | 0.00022 | 0.00011741 120 | -0.00021 | 0.000113107 121 | 0.0002 | 0.000108987 122 | -0.0002 | 0.000104995 123 | 0.00019 | 0.000101169 124 | -0.00018 | 9.74636e-05 125 | 0.00018 | 9.39119e-05 126 | -0.00017 | 9.04727e-05 127 | 0.00016 | 8.71752e-05 128 | -0.00016 | 8.39833e-05 129 | 0.00015 | 8.09218e-05 130 | -0.00015 | 7.79594e-05 131 | 0.00014 | 7.5117e-05 132 | -0.00014 | 7.23675e-05 133 | 0.00013 | 6.97286e-05 134 | -0.00013 | 6.71766e-05 135 | 0.00012 | 6.47267e-05 136 | -0.00012 | 6.23582e-05 137 | 0.00011 | 6.00837e-05 138 | -0.00011 | 5.78853e-05 139 | 0.0001 | 5.57737e-05 140 | -0.0001 | 5.37332e-05 141 | 9.7e-5 | 5.17729e-05 142 | -9.4e-5 | 4.9879e-05 143 | 9e-5 | 4.80591e-05 144 | -8.7e-5 | 4.63012e-05 145 | 8.4e-5 | 4.46117e-05 146 | -8.1e-5 | 4.298e-05 147 | 7.8e-5 | 4.14116e-05 148 | -7.5e-5 | 3.98971e-05 149 | 7.2e-5 | 3.84411e-05 150 | -7e-5 | 3.70353e-05 151 | 6.7e-5 | 3.56836e-05 152 | -6.5e-5 | 3.43788e-05 153 | 6.2e-5 | 3.31239e-05 154 | -6e-5 | 3.19128e-05 155 | 5.8e-5 | 3.07479e-05 156 | -5.6e-5 | 2.96237e-05 157 | 5.4e-5 | 2.85423e-05 158 | -5.2e-5 | 2.74988e-05 159 | 5e-5 | 2.64949e-05 160 | -4.8e-5 | 2.55263e-05 161 | 4.6e-5 | 2.45944e-05 162 | -4.4e-5 | 2.36953e-05 163 | 4.3e-5 | 2.28302e-05 164 | -4.1e-5 | 2.19956e-05 165 | 4e-5 | 2.11925e-05 166 | -3.8e-5 | 2.04179e-05 167 | 3.7e-5 | 1.96724e-05 168 | -3.6e-5 | 1.89533e-05 169 | 3.4e-5 | 1.82612e-05 170 | -3.3e-5 | 1.75938e-05 171 | 3.2e-5 | 1.69513e-05 172 | -3.1e-5 | 1.63318e-05 173 | 3e-5 | 1.57354e-05 174 | -2.8e-5 | 1.51603e-05 175 | 2.7e-5 | 1.46067e-05 176 | -2.6e-5 | 1.40728e-05 177 | 2.5e-5 | 1.35589e-05 178 | -2.5e-5 | 1.30634e-05 179 | 2.4e-5 | 1.25863e-05 180 | -2.3e-5 | 1.21263e-05 181 | 2.2e-5 | 1.16835e-05 182 | -2.1e-5 | 1.12565e-05 183 | 2e-5 | 1.08454e-05 184 | -2e-5 | 1.04491e-05 185 | 1.9e-5 | 1.00674e-05 186 | -1.8e-5 | 9.69956e-06 187 | 1.8e-5 | 9.3453e-06 188 | -1.7e-5 | 9.00381e-06 189 | 1.6e-5 | 8.67495e-06 190 | -1.6e-5 | 8.35796e-06 191 | 1.5e-5 | 8.05269e-06 192 | -1.5e-5 | 7.75844e-06 193 | 1.4e-5 | 7.47506e-06 194 | -1.4e-5 | 7.20192e-06 195 | 1.3e-5 | 6.93886e-06 196 | -1.3e-5 | 6.68532e-06 197 | 1.2e-5 | 6.44113e-06 198 | -1.2e-5 | 6.20578e-06 199 | 1.1e-5 | 5.9791e-06 200 | -1.1e-5 | 5.76063e-06 201 | 1e-5 | 5.55021e-06 202 | -1e-5 | 5.34742e-06 203 | 9.7e-6 | 5.15209e-06 204 | -9.3e-6 | 4.96385e-06 205 | 9e-6 | 4.78253e-06 206 | -8.7e-6 | 4.60779e-06 207 | 8.3e-6 | 4.43947e-06 208 | -8e-6 | 4.27727e-06 209 | 7.7e-6 | 4.12102e-06 210 | -7.5e-6 | 3.97045e-06 211 | 7.2e-6 | 3.82542e-06 212 | -6.9e-6 | 3.68565e-06 213 | 6.7e-6 | 3.55102e-06 214 | -6.4e-6 | 3.42128e-06 215 | 6.2e-6 | 3.2963e-06 216 | -6e-6 | 3.17587e-06 217 | 5.7e-6 | 3.05985e-06 218 | -5.5e-6 | 2.94806e-06 219 | 5.3e-6 | 2.84037e-06 220 | -5.1e-6 | 2.73659e-06 221 | 5e-6 | 2.63662e-06 222 | -4.8e-6 | 2.54029e-06 223 | 4.6e-6 | 2.4475e-06 224 | -4.4e-6 | 2.35808e-06 225 | 4.3e-6 | 2.27193e-06 226 | -4.1e-6 | 2.18893e-06 227 | 4e-6 | 2.10897e-06 228 | -3.8e-6 | 2.03192e-06 229 | 3.7e-6 | 1.95769e-06 230 | -3.5e-6 | 1.88617e-06 231 | 3.4e-6 | 1.81726e-06 232 | -3.3e-6 | 1.75087e-06 233 | 3.2e-6 | 1.68691e-06 234 | -3.1e-6 | 1.62528e-06 235 | 2.9e-6 | 1.5659e-06 236 | -2.8e-6 | 1.5087e-06 237 | 2.7e-6 | 1.45358e-06 238 | -2.6e-6 | 1.40048e-06 239 | 2.5e-6 | 1.34931e-06 240 | -2.4e-6 | 1.30002e-06 241 | 2.4e-6 | 1.25253e-06 242 | -2.3e-6 | 1.20677e-06 243 | 2.2e-6 | 1.16268e-06 244 | -2.1e-6 | 1.1202e-06 245 | 2e-6 | 1.07928e-06 246 | -2e-6 | 1.03985e-06 247 | 1.9e-6 | 1.00186e-06 248 | -1.8e-6 | 9.65262e-07 249 | 1.7e-6 | 9.29999e-07 250 | -1.7e-6 | 8.96023e-07 251 | 1.6e-6 | 8.63289e-07 252 | -1.6e-6 | 8.3175e-07 253 | 1.5e-6 | 8.01365e-07 254 | -1.4e-6 | 7.72088e-07 255 | 1.4e-6 | 7.43882e-07 256 | -1.3e-6 | 7.16706e-07 257 | 1.3e-6 | 6.90523e-07 258 | -1.2e-6 | 6.65296e-07 259 | 1.2e-6 | 6.40991e-07 260 | -1.2e-6 | 6.17573e-07 261 | 1.1e-6 | 5.95012e-07 262 | -1.1e-6 | 5.73274e-07 263 | 1e-6 | 5.52331e-07 264 | -1e-6 | 5.32153e-07 265 | 9.6e-7 | 5.12712e-07 266 | -9.3e-7 | 4.93981e-07 267 | 8.9e-7 | 4.75935e-07 268 | -8.6e-7 | 4.58547e-07 269 | 8.3e-7 | 4.41796e-07 270 | -8e-7 | 4.25655e-07 271 | 7.7e-7 | 4.10105e-07 272 | -7.4e-7 | 3.95123e-07 273 | 7.1e-7 | 3.80688e-07 274 | -6.9e-7 | 3.6678e-07 275 | 6.6e-7 | 3.53381e-07 276 | -6.4e-7 | 3.40471e-07 277 | 6.2e-7 | 3.28033e-07 278 | -5.9e-7 | 3.16049e-07 279 | 5.7e-7 | 3.04503e-07 280 | -5.5e-7 | 2.93378e-07 281 | 5.3e-7 | 2.8266e-07 282 | -5.1e-7 | 2.72334e-07 283 | 4.9e-7 | 2.62385e-07 284 | -4.7e-7 | 2.52799e-07 285 | 4.6e-7 | 2.43564e-07 286 | -4.4e-7 | 2.34666e-07 287 | 4.2e-7 | 2.26093e-07 288 | -4.1e-7 | 2.17833e-07 289 | 3.9e-7 | 2.09875e-07 290 | -3.8e-7 | 2.02208e-07 291 | 3.7e-7 | 1.9482e-07 292 | -3.5e-7 | 1.87703e-07 293 | 3.4e-7 | 1.80846e-07 294 | -3.3e-7 | 1.74239e-07 295 | 3.2e-7 | 1.67873e-07 296 | -3e-7 | 1.61741e-07 297 | 2.9e-7 | 1.55832e-07 298 | -2.8e-7 | 1.50139e-07 299 | 2.7e-7 | 1.44654e-07 300 | -2.6e-7 | 1.39369e-07 301 | 2.5e-7 | 1.34278e-07 302 | -2.4e-7 | 1.29372e-07 303 | 2.3e-7 | 1.24646e-07 304 | -2.3e-7 | 1.20092e-07 305 | 2.2e-7 | 1.15705e-07 306 | -2.1e-7 | 1.11478e-07 307 | 2e-7 | 1.07405e-07 308 | -1.9e-7 | 1.03481e-07 309 | 1.9e-7 | 9.9701e-08 310 | -1.8e-7 | 9.60586e-08 311 | 1.7e-7 | 9.25493e-08 312 | -1.7e-7 | 8.91683e-08 313 | 1.6e-7 | 8.59107e-08 314 | -1.6e-7 | 8.27721e-08 315 | 1.5e-7 | 7.97483e-08 316 | -1.4e-7 | 7.68348e-08 317 | 1.4e-7 | 7.40278e-08 318 | -1.3e-7 | 7.13234e-08 319 | 1.3e-7 | 6.87178e-08 320 | -1.2e-7 | 6.62073e-08 321 | 1.2e-7 | 6.37886e-08 322 | -1.2e-7 | 6.14582e-08 323 | 1.1e-7 | 5.9213e-08 324 | -1.1e-7 | 5.70497e-08 325 | 1e-7 | 5.49656e-08 326 | -9.9e-8 | 5.29575e-08 327 | 9.6e-8 | 5.10228e-08 328 | -9.2e-8 | 4.91588e-08 329 | 8.9e-8 | 4.73629e-08 330 | -8.6e-8 | 4.56326e-08 331 | 8.3e-8 | 4.39655e-08 332 | -8e-8 | 4.23594e-08 333 | 7.7e-8 | 4.08119e-08 334 | -7.4e-8 | 3.93209e-08 335 | 7.1e-8 | 3.78844e-08 336 | -6.9e-8 | 3.65004e-08 337 | 6.6e-8 | 3.51669e-08 338 | -6.4e-8 | 3.38822e-08 339 | 6.1e-8 | 3.26444e-08 340 | -5.9e-8 | 3.14518e-08 341 | 5.7e-8 | 3.03027e-08 342 | -5.5e-8 | 2.91957e-08 343 | 5.3e-8 | 2.81291e-08 344 | -5.1e-8 | 2.71015e-08 345 | 4.9e-8 | 2.61114e-08 346 | -4.7e-8 | 2.51575e-08 347 | 4.6e-8 | 2.42384e-08 348 | -4.4e-8 | 2.33529e-08 349 | 4.2e-8 | 2.24997e-08 350 | -4.1e-8 | 2.16778e-08 351 | 3.9e-8 | 2.08858e-08 352 | -3.8e-8 | 2.01228e-08 353 | 3.6e-8 | 1.93877e-08 354 | -3.5e-8 | 1.86794e-08 355 | 3.4e-8 | 1.7997e-08 356 | -3.3e-8 | 1.73395e-08 357 | 3.1e-8 | 1.6706e-08 358 | -3e-8 | 1.60957e-08 359 | 2.9e-8 | 1.55077e-08 360 | -2.8e-8 | 1.49412e-08 361 | 2.7e-8 | 1.43953e-08 362 | -2.6e-8 | 1.38694e-08 363 | 2.5e-8 | 1.33627e-08 364 | -2.4e-8 | 1.28745e-08 365 | 2.3e-8 | 1.24042e-08 366 | -2.2e-8 | 1.1951e-08 367 | 2.2e-8 | 1.15144e-08 368 | -2.1e-8 | 1.10938e-08 369 | 2e-8 | 1.06885e-08 370 | -1.9e-8 | 1.0298e-08 371 | 1.9e-8 | 9.9218e-09 372 | -1.8e-8 | 9.55933e-09 373 | 1.7e-8 | 9.2101e-09 374 | -1.7e-8 | 8.87363e-09 375 | 1.6e-8 | 8.54945e-09 376 | -1.5e-8 | 8.23712e-09 377 | 1.5e-8 | 7.93619e-09 378 | -1.4e-8 | 7.64626e-09 379 | 1.4e-8 | 7.36692e-09 380 | -1.3e-8 | 7.09779e-09 381 | 1.3e-8 | 6.83849e-09 382 | -1.2e-8 | 6.58866e-09 383 | 1.2e-8 | 6.34796e-09 384 | -1.1e-8 | 6.11605e-09 385 | 1.1e-8 | 5.89261e-09 386 | -1.1e-8 | 5.67734e-09 387 | 1e-8 | 5.46993e-09 388 | -9.9e-9 | 5.2701e-09 389 | 9.5e-9 | 5.07757e-09 390 | -9.2e-9 | 4.89207e-09 391 | 8.9e-9 | 4.71335e-09 392 | -8.5e-9 | 4.54116e-09 393 | 8.2e-9 | 4.37526e-09 394 | -7.9e-9 | 4.21541e-09 395 | 7.6e-9 | 4.06141e-09 396 | -7.3e-9 | 3.91304e-09 397 | 7.1e-9 | 3.77008e-09 398 | -6.8e-9 | 3.63235e-09 399 | 6.6e-9 | 3.49965e-09 400 | -6.3e-9 | 3.3718e-09 401 | 6.1e-9 | 3.24862e-09 402 | -5.9e-9 | 3.12994e-09 403 | 5.7e-9 | 3.01559e-09 404 | -5.5e-9 | 2.90543e-09 405 | 5.3e-9 | 2.79928e-09 406 | -5.1e-9 | 2.69702e-09 407 | 4.9e-9 | 2.59849e-09 408 | -4.7e-9 | 2.50356e-09 409 | 4.5e-9 | 2.4121e-09 410 | -4.4e-9 | 2.32398e-09 411 | 4.2e-9 | 2.23907e-09 412 | -4.1e-9 | 2.15727e-09 413 | 3.9e-9 | 2.07846e-09 414 | -3.8e-9 | 2.00253e-09 415 | 3.6e-9 | 1.92937e-09 416 | -3.5e-9 | 1.85889e-09 417 | 3.4e-9 | 1.79098e-09 418 | -3.2e-9 | 1.72555e-09 419 | 3.1e-9 | 1.66251e-09 420 | -3e-9 | 1.60177e-09 421 | 2.9e-9 | 1.54326e-09 422 | -2.8e-9 | 1.48688e-09 423 | 2.7e-9 | 1.43256e-09 424 | -2.6e-9 | 1.38022e-09 425 | 2.5e-9 | 1.3298e-09 426 | -2.4e-9 | 1.28122e-09 427 | 2.3e-9 | 1.23441e-09 428 | -2.2e-9 | 1.18931e-09 429 | 2.2e-9 | 1.14587e-09 430 | -2.1e-9 | 1.104e-09 431 | 2e-9 | 1.06367e-09 432 | -1.9e-9 | 1.02481e-09 433 | 1.9e-9 | 9.87373e-10 434 | -1.8e-9 | 9.51302e-10 435 | 1.7e-9 | 9.16548e-10 436 | -1.7e-9 | 8.83064e-10 437 | 1.6e-9 | 8.50803e-10 438 | -1.5e-9 | 8.19721e-10 439 | 1.5e-9 | 7.89774e-10 440 | -1.4e-9 | 7.60921e-10 441 | 1.4e-9 | 7.33123e-10 442 | -1.3e-9 | 7.0634e-10 443 | 1.3e-9 | 6.80535e-10 444 | -1.2e-9 | 6.55673e-10 445 | 1.2e-9 | 6.3172e-10 446 | -1.1e-9 | 6.08641e-10 447 | 1.1e-9 | 5.86406e-10 448 | -1.1e-9 | 5.64983e-10 449 | 1e-9 | 5.44343e-10 450 | -9.8e-10 | 5.24456e-10 451 | 9.5e-10 | 5.05296e-10 452 | -9.1e-10 | 4.86837e-10 453 | 8.8e-10 | 4.69051e-10 454 | -8.5e-10 | 4.51915e-10 455 | 8.2e-10 | 4.35405e-10 456 | -7.9e-10 | 4.19499e-10 457 | 7.6e-10 | 4.04174e-10 458 | -7.3e-10 | 3.89408e-10 459 | 7e-10 | 3.75182e-10 460 | -6.8e-10 | 3.61475e-10 461 | 6.5e-10 | 3.4827e-10 462 | -6.3e-10 | 3.35547e-10 463 | 6.1e-10 | 3.23288e-10 464 | -5.8e-10 | 3.11477e-10 465 | 5.6e-10 | 3.00098e-10 466 | -5.4e-10 | 2.89134e-10 467 | 5.2e-10 | 2.78572e-10 468 | -5e-10 | 2.68395e-10 469 | 4.9e-10 | 2.5859e-10 470 | -4.7e-10 | 2.49143e-10 471 | 4.5e-10 | 2.40041e-10 472 | -4.3e-10 | 2.31271e-10 473 | 4.2e-10 | 2.22822e-10 474 | -4e-10 | 2.14682e-10 475 | 3.9e-10 | 2.06839e-10 476 | -3.7e-10 | 1.99283e-10 477 | 3.6e-10 | 1.92002e-10 478 | -3.5e-10 | 1.84988e-10 479 | 3.3e-10 | 1.7823e-10 480 | -3.2e-10 | 1.71719e-10 481 | 3.1e-10 | 1.65445e-10 482 | -3e-10 | 1.59401e-10 483 | 2.9e-10 | 1.53578e-10 484 | -2.8e-10 | 1.47967e-10 485 | 2.7e-10 | 1.42562e-10 486 | -2.6e-10 | 1.37354e-10 487 | 2.5e-10 | 1.32336e-10 488 | -2.4e-10 | 1.27501e-10 489 | 2.3e-10 | 1.22843e-10 490 | -2.2e-10 | 1.18355e-10 491 | 2.1e-10 | 1.14031e-10 492 | -2.1e-10 | 1.09866e-10 493 | 2e-10 | 1.05852e-10 494 | -1.9e-10 | 1.01985e-10 495 | 1.8e-10 | 9.82592e-11 > ## slow convergence > stopifnot(okFit(rm2)) > > > ###>> 2 ------------------- family = binomial ----------------------------------- > > ## Using *factor* y ... > x <- seq(0,5, length = 120) > summary(px <- plogis(-5 + 2*x)) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.006693 0.075897 0.500000 0.500000 0.924103 0.993307 > set.seed(7) > (f <- factor(rbinom(length(x), 1, prob=px))) [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0 [38] 0 0 0 0 0 0 0 1 0 0 1 1 1 0 1 1 0 1 0 0 0 0 0 1 1 0 0 1 0 1 1 1 1 1 1 0 1 [75] 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 [112] 1 1 1 1 1 1 1 1 1 Levels: 0 1 > > summary(m.c0 <- glm (f ~ x, family = binomial)) Call: glm(formula = f ~ x, family = binomial) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -3.9428 0.7259 -5.431 5.59e-08 *** x 1.6943 0.2869 5.906 3.50e-09 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 165.822 on 119 degrees of freedom Residual deviance: 87.395 on 118 degrees of freedom AIC: 91.395 Number of Fisher Scoring iterations: 5 > summary(m.r0 <- glmrob(f ~ x, family = binomial)) Call: glmrob(formula = f ~ x, family = binomial) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -3.8800 0.7779 -4.988 6.10e-07 *** x 1.6341 0.3044 5.369 7.91e-08 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Robustness weights w.r * w.x: 110 weights are ~= 1. The remaining 10 ones are 19 23 25 29 45 48 73 80 81 82 0.3586 0.4113 0.4406 0.5054 0.8754 0.9704 0.7903 0.6215 0.6005 0.5802 Number of observations: 120 Fitted by method 'Mqle' (in 4 iterations) (Dispersion parameter for binomial family taken to be 1) No deviance values available Algorithmic parameters: acc tcc 0.0001 1.3450 maxit 50 test.acc "coef" > > ## add outliers --- in y: > f. <- f > f.[i1 <- 2:3] <- 1 > f.[i0 <- 110+c(1,7)] <- 0 > m.c1 <- glm (f. ~ x, family = binomial) > summary(m.r1 <- glmrob(f. ~ x, family = binomial)) ## hmm, not so robust? Call: glmrob(formula = f. ~ x, family = binomial) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -3.1082 0.6269 -4.958 7.12e-07 *** x 1.3162 0.2402 5.479 4.29e-08 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Robustness weights w.r * w.x: 107 weights are ~= 1. The remaining 13 ones are summarized as Min. 1st Qu. Median Mean 3rd Qu. Max. 0.2574 0.3039 0.5521 0.5563 0.6966 0.9598 Number of observations: 120 Fitted by method 'Mqle' (in 5 iterations) (Dispersion parameter for binomial family taken to be 1) No deviance values available Algorithmic parameters: acc tcc 0.0001 1.3450 maxit 50 test.acc "coef" > stopifnot(m.r1$w.r[c(i0,i1)] < 1/3, # well, at least down weighted + ## and coefficients change less : + (coef(m.r1) - coef(m.c0)) / (coef(m.c1) - coef(m.c0)) < 1) > assert.EQ(c("(Intercept)" = -3.10817337603974, x = 1.31618564057790), + coef(m.r1), tol= 1e-14, giveRE=TRUE) Mean relative difference: 1.687031e-15 > > y <- as.numeric(as.character(f.)) > m.r2 <- BYlogreg(x0=x, y=y, trace=TRUE, maxhalf= 10) k= 1, s1= 0.34255783: => new s1= 0.31673535 k= 2, s1= 0.31673535: => new s1= 0.30965429 k= 3, s1= 0.30965429: => new s1= 0.31134409 Convergence Achieved > m.r2A <- BYlogreg(x0=x, y=y, trace= 2 , maxhalf= 15) k= 1, s1= 0.34255783: => new s1= 0.31673535, obj= 0.544734224: jh= 1, hstep= 0.5 => new obj= 0.968449422 jh= 2, hstep= 0.25 => new obj= 0.721112749 jh= 3, hstep= 0.125 => new obj= 0.582862929 jh= 4, hstep= 0.0625 => new obj= 0.552239517 jh= 5, hstep= 0.03125 => new obj= 0.546125983 jh= 6, hstep= 0.015625 => new obj= 0.544900448 jh= 7, hstep= 0.0078125 => new obj= 0.544692433 k= 2, s1= 0.31673535: => new s1= 0.30965429, obj= 0.544674745: jh= 1, hstep= 0.5 => new obj= 0.699044068 jh= 2, hstep= 0.25 => new obj= 0.611232275 jh= 3, hstep= 0.125 => new obj= 0.567137872 jh= 4, hstep= 0.0625 => new obj= 0.551165685 jh= 5, hstep= 0.03125 => new obj= 0.546354383 jh= 6, hstep= 0.015625 => new obj= 0.545068151 jh= 7, hstep= 0.0078125 => new obj= 0.544752043 jh= 8, hstep= 0.00390625 => new obj= 0.54468258 jh= 9, hstep= 0.001953125 => new obj= 0.544670853 k= 3, s1= 0.30965429: => new s1= 0.31134409, obj= 0.544669825: jh= 1, hstep= 0.5 => new obj= 0.977985364 jh= 2, hstep= 0.25 => new obj= 0.732177489 jh= 3, hstep= 0.125 => new obj= 0.587549933 jh= 4, hstep= 0.0625 => new obj= 0.553913233 jh= 5, hstep= 0.03125 => new obj= 0.5467962 jh= 6, hstep= 0.015625 => new obj= 0.545173095 jh= 7, hstep= 0.0078125 => new obj= 0.544789191 jh= 8, hstep= 0.00390625 => new obj= 0.544697478 jh= 9, hstep= 0.001953125 => new obj= 0.544675771 jh=10, hstep= 0.0009765625 => new obj= 0.544670839 jh=11, hstep=0.00048828125 => new obj= 0.544669844 jh=12, hstep=0.00024414062 => new obj= 0.544669713 k= 4, s1= 0.31134409: => new s1= 0.31112977, obj= 0.544669696: jh= 1, hstep= 0.5 => new obj= 0.978373379 jh= 2, hstep= 0.25 => new obj= 0.732641189 jh= 3, hstep= 0.125 => new obj= 0.58775415 jh= 4, hstep= 0.0625 => new obj= 0.553988479 jh= 5, hstep= 0.03125 => new obj= 0.546828008 jh= 6, hstep= 0.015625 => new obj= 0.545187648 jh= 7, hstep= 0.0078125 => new obj= 0.544796093 jh= 8, hstep= 0.00390625 => new obj= 0.544700782 jh= 9, hstep= 0.001953125 => new obj= 0.544677343 jh=10, hstep= 0.0009765625 => new obj= 0.544671557 jh=11, hstep=0.00048828125 => new obj= 0.544670138 jh=12, hstep=0.00024414062 => new obj= 0.544669795 jh=13, hstep=0.00012207031 => new obj= 0.544669715 jh=14, hstep=6.1035156e-05 => new obj= 0.544669698 jh=15, hstep=3.0517578e-05 => new obj= 0.544669695 k= 5, s1= 0.31112977: => new s1= 0.31110309, obj= 0.544669695: jh= 1, hstep= 0.5 => new obj= 0.699105497 jh= 2, hstep= 0.25 => new obj= 0.611555159 jh= 3, hstep= 0.125 => new obj= 0.567491276 jh= 4, hstep= 0.0625 => new obj= 0.551421622 jh= 5, hstep= 0.03125 => new obj= 0.546506517 jh= 6, hstep= 0.015625 => new obj= 0.545149258 jh= 7, hstep= 0.0078125 => new obj= 0.544792088 jh= 8, hstep= 0.00390625 => new obj= 0.544700561 jh= 9, hstep= 0.001953125 => new obj= 0.54467743 jh=10, hstep= 0.0009765625 => new obj= 0.544671626 jh=11, hstep=0.00048828125 => new obj= 0.544670175 jh=12, hstep=0.00024414062 => new obj= 0.544669814 jh=13, hstep=0.00012207031 => new obj= 0.544669724 jh=14, hstep=6.1035156e-05 => new obj= 0.544669702 jh=15, hstep=3.0517578e-05 => new obj= 0.544669697 Convergence Achieved > ## different.. but not so much: > iB <- 1:5 > assert.EQ(m.r2A[iB], m.r2[iB], tol = .003, giveRE=TRUE) Component "objective": Mean relative difference: 2.378848e-07 Component "coefficients": Mean relative difference: 0.0006633958 Component "cov": Mean relative difference: 0.002200778 Component "sterror": Mean relative difference: 0.001061053 > > > assert.EQ(c("Intercept" = -2.9554950286, x = 1.2574679132), + ## 32-bit{ada-5} -2.95549502890363 1.25746791332613 + m.r2$coef, tol=8e-10, giveRE=TRUE)# seen 5.316e-10 for --disable-long-double Mean relative difference: 5.061879e-10 > assert.EQ( c(0.685919891749065, 0.256419206157062), + ## 32-bit{ada-5}: + ## 0.685919891858219, 0.256419206203016) + m.r2$sterror, tol=4e-9)# seen 1.025e-9 for --disable-long-double > > data(foodstamp) > str(foodstamp) 'data.frame': 150 obs. of 4 variables: $ participation: int 0 0 0 0 0 0 0 0 0 0 ... $ tenancy : int 1 1 1 1 0 1 1 1 0 1 ... $ suppl.income : int 0 0 1 0 0 0 0 0 0 0 ... $ income : int 271 287 714 521 0 518 458 1266 350 168 ... > ## Model with 'income' instead of log(income+1) is "interesting" > ## because BYlogreg() needs maxhalf > 10 for convergence! > m.fs0 <- glm (participation ~ ., family=binomial, data=foodstamp) > m.fs0QL <- glmrob(participation ~ ., family=binomial, data=foodstamp) > y.fs <- foodstamp[,"participation"] > X.fs0 <- model.matrix(m.fs0) > head(X.fs0) (Intercept) tenancy suppl.income income 1 1 1 0 271 2 1 1 0 287 3 1 1 1 714 4 1 1 0 521 5 1 0 0 0 6 1 1 0 518 > ## (former default) maxhalf = 10 leads to too early convergence: > m.fsWBY. <- BYlogreg(x0=X.fs0, y=y.fs, + addIntercept=FALSE, trace=TRUE, maxhalf=10) k= 1, s1= 0.51017267: => new s1= 0.56467262 Convergence Achieved > m.fs.BY. <- BYlogreg(x0=X.fs0, y=y.fs, initwml=FALSE, + addIntercept=FALSE, trace=TRUE, maxhalf=10) k= 1, s1= 0.5119608: => new s1= 0.51098755 Convergence Achieved > m.fsWBY <- BYlogreg(x0=X.fs0, y=y.fs, + addIntercept=FALSE, trace=TRUE, maxhalf=18) k= 1, s1= 0.51017267: => new s1= 0.56467262 k= 2, s1= 0.56467262: => new s1= 0.55168008 k= 3, s1= 0.55168008: => new s1= 0.54529145 k= 4, s1= 0.54529145: => new s1= 0.54370564 k= 5, s1= 0.54370564: => new s1= 0.54291823 Convergence Achieved > m.fs.BY <- BYlogreg(x0=X.fs0, y=y.fs, initwml=FALSE, + addIntercept=FALSE, trace=TRUE, maxhalf=18) k= 1, s1= 0.5119608: => new s1= 0.51098755 k= 2, s1= 0.51098755: => new s1= 0.52513485 k= 3, s1= 0.52513485: => new s1= 0.53948806 k= 4, s1= 0.53948806: => new s1= 0.54129516 k= 5, s1= 0.54129516: => new s1= 0.54310324 Convergence Achieved > > assert.EQ(m.fsWBY.[iB], m.fsWBY[iB], tol= 0.07)## almost 7% different > assert.EQ(m.fs.BY.[iB], m.fs.BY[iB], tol= 0.08) > > foodSt <- within(foodstamp, { logInc <- log(1 + income) ; rm(income) }) > > m.fsML <- glm (participation ~ ., family=binomial, data=foodSt) > m.fsQL <- glmrob(participation ~ ., family=binomial, data=foodSt) > X.fs <- model.matrix(m.fsML) > stopifnot(dim(X.fs) == c(150, 4)) # including intercept! > try(## FIXME -- Mahalanobis fails with singular matrix, here: + m.fsWBY <- BYlogreg(x0=X.fs, y=y.fs, + addIntercept=FALSE, trace=TRUE, maxhalf=18) + ) Error in .MCDsingularityMsg(ans$singularity, ans$n.obs) : illegal 'singularity$kind' > ## maxhalf=18 is too much --> no convergence (in 1000 steps) > m.fs.BY <- BYlogreg(x0=X.fs, y=y.fs, initwml=FALSE, + addIntercept=FALSE, trace=TRUE, maxhalf=18) k= 1, s1= 0.43873492: => new s1= 0.45207408 k= 2, s1= 0.45207408: => new s1= 0.45506753 k= 3, s1= 0.45506753: => new s1= 0.45799437 k= 4, s1= 0.45799437: => new s1= 0.46031985 k= 5, s1= 0.46031985: => new s1= 0.45881914 k= 6, s1= 0.45881914: => new s1= 0.45751194 k= 7, s1= 0.45751194: => new s1= 0.45932753 k= 8, s1= 0.45932753: => new s1= 0.45785407 k= 9, s1= 0.45785407: => new s1= 0.45710539 k=10, s1= 0.45710539: => new s1= 0.45793621 k=11, s1= 0.45793621: => new s1= 0.45742459 k=12, s1= 0.45742459: => new s1= 0.45768312 k=13, s1= 0.45768312: => new s1= 0.45691269 k=14, s1= 0.45691269: => new s1= 0.4561501 k=15, s1= 0.4561501: => new s1= 0.455997 k=16, s1= 0.455997: => new s1= 0.4562702 k=17, s1= 0.4562702: => new s1= 0.45665261 k=18, s1= 0.45665261: => new s1= 0.45672591 k=19, s1= 0.45672591: => new s1= 0.45711865 k=20, s1= 0.45711865: => new s1= 0.45692986 k=21, s1= 0.45692986: => new s1= 0.4568429 k=22, s1= 0.4568429: => new s1= 0.45678798 k=23, s1= 0.45678798: => new s1= 0.45690093 k=24, s1= 0.45690093: => new s1= 0.45682335 k=25, s1= 0.45682335: => new s1= 0.45663851 k=26, s1= 0.45663851: => new s1= 0.45702503 k=27, s1= 0.45702503: => new s1= 0.45721719 k=28, s1= 0.45721719: => new s1= 0.45730701 k=29, s1= 0.45730701: => new s1= 0.45740594 k=30, s1= 0.45740594: => new s1= 0.45733832 k=31, s1= 0.45733832: => new s1= 0.45744587 k=32, s1= 0.45744587: => new s1= 0.45737153 k=33, s1= 0.45737153: => new s1= 0.45766256 k=34, s1= 0.45766256: => new s1= 0.45757041 k=35, s1= 0.45757041: => new s1= 0.45748704 k=36, s1= 0.45748704: => new s1= 0.45724295 k=37, s1= 0.45724295: => new s1= 0.45743597 k=38, s1= 0.45743597: => new s1= 0.45762776 k=39, s1= 0.45762776: => new s1= 0.45770141 k=40, s1= 0.45770141: => new s1= 0.45765689 k=41, s1= 0.45765689: => new s1= 0.45770888 k=42, s1= 0.45770888: => new s1= 0.45762472 k=43, s1= 0.45762472: => new s1= 0.45766739 k=44, s1= 0.45766739: => new s1= 0.45772587 k=45, s1= 0.45772587: => new s1= 0.45768048 k=46, s1= 0.45768048: => new s1= 0.45773336 k=47, s1= 0.45773336: => new s1= 0.45764719 k=48, s1= 0.45764719: => new s1= 0.4576901 k=49, s1= 0.4576901: => new s1= 0.45775023 k=50, s1= 0.45775023: => new s1= 0.45770378 k=51, s1= 0.45770378: => new s1= 0.45775757 k=52, s1= 0.45775757: => new s1= 0.45771345 k=53, s1= 0.45771345: => new s1= 0.45776641 k=54, s1= 0.45776641: => new s1= 0.45772245 k=55, s1= 0.45772245: => new s1= 0.4577756 k=56, s1= 0.4577756: => new s1= 0.45773141 k=57, s1= 0.45773141: => new s1= 0.45778476 k=58, s1= 0.45778476: => new s1= 0.45774021 k=59, s1= 0.45774021: => new s1= 0.45779394 k=60, s1= 0.45779394: => new s1= 0.45774897 k=61, s1= 0.45774897: => new s1= 0.45780308 k=62, s1= 0.45780308: => new s1= 0.45775772 k=63, s1= 0.45775772: => new s1= 0.45781219 k=64, s1= 0.45781219: => new s1= 0.45776641 k=65, s1= 0.45776641: => new s1= 0.45782132 k=66, s1= 0.45782132: => new s1= 0.45777513 k=67, s1= 0.45777513: => new s1= 0.45783033 k=68, s1= 0.45783033: => new s1= 0.45778372 k=69, s1= 0.45778372: => new s1= 0.4578393 k=70, s1= 0.4578393: => new s1= 0.45779236 k=71, s1= 0.45779236: => new s1= 0.45784827 k=72, s1= 0.45784827: => new s1= 0.45780096 k=73, s1= 0.45780096: => new s1= 0.45785723 k=74, s1= 0.45785723: => new s1= 0.45780953 k=75, s1= 0.45780953: => new s1= 0.45786609 k=76, s1= 0.45786609: => new s1= 0.45781804 k=77, s1= 0.45781804: => new s1= 0.45787496 k=78, s1= 0.45787496: => new s1= 0.45782653 k=79, s1= 0.45782653: => new s1= 0.4578838 k=80, s1= 0.4578838: => new s1= 0.45783503 k=81, s1= 0.45783503: => new s1= 0.45789261 k=82, s1= 0.45789261: => new s1= 0.45784346 k=83, s1= 0.45784346: => new s1= 0.45790136 k=84, s1= 0.45790136: => new s1= 0.45785188 k=85, s1= 0.45785188: => new s1= 0.45791012 k=86, s1= 0.45791012: => new s1= 0.4578603 k=87, s1= 0.4578603: => new s1= 0.45791879 k=88, s1= 0.45791879: => new s1= 0.45786866 k=89, s1= 0.45786866: => new s1= 0.4579275 k=90, s1= 0.4579275: => new s1= 0.45787698 k=91, s1= 0.45787698: => new s1= 0.45793614 k=92, s1= 0.45793614: => new s1= 0.45788528 k=93, s1= 0.45788528: => new s1= 0.45794472 k=94, s1= 0.45794472: => new s1= 0.45789354 k=95, s1= 0.45789354: => new s1= 0.45795333 k=96, s1= 0.45795333: => new s1= 0.45790181 k=97, s1= 0.45790181: => new s1= 0.45796184 k=98, s1= 0.45796184: => new s1= 0.45791007 k=99, s1= 0.45791007: => new s1= 0.45797039 k=100, s1= 0.45797039: => new s1= 0.45791826 k=101, s1= 0.45791826: => new s1= 0.45797891 k=102, s1= 0.45797891: => new s1= 0.45792648 k=103, s1= 0.45792648: => new s1= 0.45798736 k=104, s1= 0.45798736: => new s1= 0.45793457 k=105, s1= 0.45793457: => new s1= 0.45799581 k=106, s1= 0.45799581: => new s1= 0.45794274 k=107, s1= 0.45794274: => new s1= 0.45800421 k=108, s1= 0.45800421: => new s1= 0.45795086 k=109, s1= 0.45795086: => new s1= 0.45801262 k=110, s1= 0.45801262: => new s1= 0.45795894 k=111, s1= 0.45795894: => new s1= 0.45802095 k=112, s1= 0.45802095: => new s1= 0.45796699 k=113, s1= 0.45796699: => new s1= 0.45802927 k=114, s1= 0.45802927: => new s1= 0.45797505 k=115, s1= 0.45797505: => new s1= 0.45803756 k=116, s1= 0.45803756: => new s1= 0.45798302 k=117, s1= 0.45798302: => new s1= 0.45804586 k=118, s1= 0.45804586: => new s1= 0.45799101 k=119, s1= 0.45799101: => new s1= 0.45805409 k=120, s1= 0.45805409: => new s1= 0.45799898 k=121, s1= 0.45799898: => new s1= 0.45806231 k=122, s1= 0.45806231: => new s1= 0.45800693 k=123, s1= 0.45800693: => new s1= 0.45807054 k=124, s1= 0.45807054: => new s1= 0.45801487 k=125, s1= 0.45801487: => new s1= 0.45807868 k=126, s1= 0.45807868: => new s1= 0.45802274 k=127, s1= 0.45802274: => new s1= 0.45808685 k=128, s1= 0.45808685: => new s1= 0.4580306 k=129, s1= 0.4580306: => new s1= 0.45809494 k=130, s1= 0.45809494: => new s1= 0.45803846 k=131, s1= 0.45803846: => new s1= 0.45810305 k=132, s1= 0.45810305: => new s1= 0.45804629 k=133, s1= 0.45804629: => new s1= 0.45811113 k=134, s1= 0.45811113: => new s1= 0.45805407 k=135, s1= 0.45805407: => new s1= 0.45811915 k=136, s1= 0.45811915: => new s1= 0.45806185 k=137, s1= 0.45806185: => new s1= 0.45812719 k=138, s1= 0.45812719: => new s1= 0.45806964 k=139, s1= 0.45806964: => new s1= 0.45813518 k=140, s1= 0.45813518: => new s1= 0.45807733 k=141, s1= 0.45807733: => new s1= 0.45814315 k=142, s1= 0.45814315: => new s1= 0.45808508 k=143, s1= 0.45808508: => new s1= 0.45811807 k=144, s1= 0.45811807: => new s1= 0.45875099 k=145, s1= 0.45875099: => new s1= 0.4588043 k=146, s1= 0.4588043: => new s1= 0.45876966 k=147, s1= 0.45876966: => new s1= 0.45891749 k=148, s1= 0.45891749: => new s1= 0.45882548 k=149, s1= 0.45882548: => new s1= 0.45876335 k=150, s1= 0.45876335: => new s1= 0.45880829 k=151, s1= 0.45880829: => new s1= 0.45886206 k=152, s1= 0.45886206: => new s1= 0.45882662 k=153, s1= 0.45882662: => new s1= 0.45897859 k=154, s1= 0.45897859: => new s1= 0.45888612 k=155, s1= 0.45888612: => new s1= 0.45884446 k=156, s1= 0.45884446: => new s1= 0.45888274 k=157, s1= 0.45888274: => new s1= 0.45877446 k=158, s1= 0.45877446: => new s1= 0.45887029 k=159, s1= 0.45887029: => new s1= 0.45891371 k=160, s1= 0.45891371: => new s1= 0.45900933 k=161, s1= 0.45900933: => new s1= 0.45896491 k=162, s1= 0.45896491: => new s1= 0.45893301 k=163, s1= 0.45893301: => new s1= 0.45900098 k=164, s1= 0.45900098: => new s1= 0.45895926 k=165, s1= 0.45895926: => new s1= 0.45887462 k=166, s1= 0.45887462: => new s1= 0.4589219 k=167, s1= 0.4589219: => new s1= 0.458966 k=168, s1= 0.458966: => new s1= 0.45902706 k=169, s1= 0.45902706: => new s1= 0.45898649 k=170, s1= 0.45898649: => new s1= 0.45885032 k=171, s1= 0.45885032: => new s1= 0.45904351 k=172, s1= 0.45904351: => new s1= 0.45909109 k=173, s1= 0.45909109: => new s1= 0.45913638 k=174, s1= 0.45913638: => new s1= 0.45917735 k=175, s1= 0.45917735: => new s1= 0.45914247 k=176, s1= 0.45914247: => new s1= 0.45918407 k=177, s1= 0.45918407: => new s1= 0.45914886 k=178, s1= 0.45914886: => new s1= 0.45919012 k=179, s1= 0.45919012: => new s1= 0.45915497 k=180, s1= 0.45915497: => new s1= 0.4591968 k=181, s1= 0.4591968: => new s1= 0.45916133 k=182, s1= 0.45916133: => new s1= 0.45920283 k=183, s1= 0.45920283: => new s1= 0.45916742 k=184, s1= 0.45916742: => new s1= 0.45920947 k=185, s1= 0.45920947: => new s1= 0.45917375 k=186, s1= 0.45917375: => new s1= 0.45921549 k=187, s1= 0.45921549: => new s1= 0.45917982 k=188, s1= 0.45917982: => new s1= 0.45922214 k=189, s1= 0.45922214: => new s1= 0.45918615 k=190, s1= 0.45918615: => new s1= 0.45922808 k=191, s1= 0.45922808: => new s1= 0.45919217 k=192, s1= 0.45919217: => new s1= 0.45923475 k=193, s1= 0.45923475: => new s1= 0.4591985 k=194, s1= 0.4591985: => new s1= 0.45924058 k=195, s1= 0.45924058: => new s1= 0.45920445 k=196, s1= 0.45920445: => new s1= 0.45924739 k=197, s1= 0.45924739: => new s1= 0.45921086 k=198, s1= 0.45921086: => new s1= 0.45925308 k=199, s1= 0.45925308: => new s1= 0.45921674 k=200, s1= 0.45921674: => new s1= 0.45925996 k=201, s1= 0.45925996: => new s1= 0.45922315 k=202, s1= 0.45922315: => new s1= 0.4592654 k=203, s1= 0.4592654: => new s1= 0.45922891 k=204, s1= 0.45922891: => new s1= 0.4592726 k=205, s1= 0.4592726: => new s1= 0.45923547 k=206, s1= 0.45923547: => new s1= 0.45927759 k=207, s1= 0.45927759: => new s1= 0.45924099 k=208, s1= 0.45924099: => new s1= 0.45928526 k=209, s1= 0.45928526: => new s1= 0.45924776 k=210, s1= 0.45924776: => new s1= 0.45928954 k=211, s1= 0.45928954: => new s1= 0.45925294 k=212, s1= 0.45925294: => new s1= 0.45929817 k=213, s1= 0.45929817: => new s1= 0.4592602 k=214, s1= 0.4592602: => new s1= 0.45930111 k=215, s1= 0.45930111: => new s1= 0.4592647 k=216, s1= 0.4592647: => new s1= 0.45931141 k=217, s1= 0.45931141: => new s1= 0.45927285 k=218, s1= 0.45927285: => new s1= 0.45931203 k=219, s1= 0.45931203: => new s1= 0.45927619 k=220, s1= 0.45927619: => new s1= 0.45928475 k=221, s1= 0.45928475: => new s1= 0.45927665 k=222, s1= 0.45927665: => new s1= 0.45930646 k=223, s1= 0.45930646: => new s1= 0.45929767 k=224, s1= 0.45929767: => new s1= 0.45930647 k=225, s1= 0.45930647: => new s1= 0.45929694 k=226, s1= 0.45929694: => new s1= 0.45930587 k=227, s1= 0.45930587: => new s1= 0.45929657 k=228, s1= 0.45929657: => new s1= 0.45930578 k=229, s1= 0.45930578: => new s1= 0.45929673 k=230, s1= 0.45929673: => new s1= 0.45930618 k=231, s1= 0.45930618: => new s1= 0.45929733 k=232, s1= 0.45929733: => new s1= 0.45930697 k=233, s1= 0.45930697: => new s1= 0.45929827 k=234, s1= 0.45929827: => new s1= 0.45930806 k=235, s1= 0.45930806: => new s1= 0.45929949 k=236, s1= 0.45929949: => new s1= 0.45930939 k=237, s1= 0.45930939: => new s1= 0.45930091 k=238, s1= 0.45930091: => new s1= 0.4593109 k=239, s1= 0.4593109: => new s1= 0.45930249 k=240, s1= 0.45930249: => new s1= 0.45931255 k=241, s1= 0.45931255: => new s1= 0.45930419 k=242, s1= 0.45930419: => new s1= 0.45933311 k=243, s1= 0.45933311: => new s1= 0.45930361 k=244, s1= 0.45930361: => new s1= 0.45931342 k=245, s1= 0.45931342: => new s1= 0.45930546 k=246, s1= 0.45930546: => new s1= 0.45933568 k=247, s1= 0.45933568: => new s1= 0.45932686 k=248, s1= 0.45932686: => new s1= 0.45933572 k=249, s1= 0.45933572: => new s1= 0.4593262 k=250, s1= 0.4593262: => new s1= 0.45933519 k=251, s1= 0.45933519: => new s1= 0.45932585 k=252, s1= 0.45932585: => new s1= 0.45933511 k=253, s1= 0.45933511: => new s1= 0.45932601 k=254, s1= 0.45932601: => new s1= 0.45933549 k=255, s1= 0.45933549: => new s1= 0.45932659 k=256, s1= 0.45932659: => new s1= 0.45933625 k=257, s1= 0.45933625: => new s1= 0.45932751 k=258, s1= 0.45932751: => new s1= 0.45933732 k=259, s1= 0.45933732: => new s1= 0.45932868 k=260, s1= 0.45932868: => new s1= 0.4593386 k=261, s1= 0.4593386: => new s1= 0.45933006 k=262, s1= 0.45933006: => new s1= 0.45934007 k=263, s1= 0.45934007: => new s1= 0.4593316 k=264, s1= 0.4593316: => new s1= 0.45936041 k=265, s1= 0.45936041: => new s1= 0.45933012 k=266, s1= 0.45933012: => new s1= 0.45934019 k=267, s1= 0.45934019: => new s1= 0.45933221 k=268, s1= 0.45933221: => new s1= 0.45936252 k=269, s1= 0.45936252: => new s1= 0.45935353 k=270, s1= 0.45935353: => new s1= 0.45936244 k=271, s1= 0.45936244: => new s1= 0.45935289 k=272, s1= 0.45935289: => new s1= 0.45936191 k=273, s1= 0.45936191: => new s1= 0.45935254 k=274, s1= 0.45935254: => new s1= 0.45936182 k=275, s1= 0.45936182: => new s1= 0.45935268 k=276, s1= 0.45935268: => new s1= 0.45936219 k=277, s1= 0.45936219: => new s1= 0.45935323 k=278, s1= 0.45935323: => new s1= 0.45936292 k=279, s1= 0.45936292: => new s1= 0.45935412 k=280, s1= 0.45935412: => new s1= 0.45936394 k=281, s1= 0.45936394: => new s1= 0.45935524 k=282, s1= 0.45935524: => new s1= 0.45936518 k=283, s1= 0.45936518: => new s1= 0.45935658 k=284, s1= 0.45935658: => new s1= 0.45936661 k=285, s1= 0.45936661: => new s1= 0.45935808 k=286, s1= 0.45935808: => new s1= 0.45938693 k=287, s1= 0.45938693: => new s1= 0.45935619 k=288, s1= 0.45935619: => new s1= 0.45938568 k=289, s1= 0.45938568: => new s1= 0.45937707 k=290, s1= 0.45937707: => new s1= 0.45937883 k=291, s1= 0.45937883: => new s1= 0.45937246 k=292, s1= 0.45937246: => new s1= 0.45937558 k=293, s1= 0.45937558: => new s1= 0.45937024 k=294, s1= 0.45937024: => new s1= 0.4593741 k=295, s1= 0.4593741: => new s1= 0.45934332 k=296, s1= 0.45934332: => new s1= 0.45938791 k=297, s1= 0.45938791: => new s1= 0.45928512 k=298, s1= 0.45928512: => new s1= 0.45947836 k=299, s1= 0.45947836: => new s1= 0.45952563 k=300, s1= 0.45952563: => new s1= 0.4595323 k=301, s1= 0.4595323: => new s1= 0.45952701 k=302, s1= 0.45952701: => new s1= 0.45956268 k=303, s1= 0.45956268: => new s1= 0.45955362 k=304, s1= 0.45955362: => new s1= 0.45955543 k=305, s1= 0.45955543: => new s1= 0.45954878 k=306, s1= 0.45954878: => new s1= 0.45955206 k=307, s1= 0.45955206: => new s1= 0.4595462 k=308, s1= 0.4595462: => new s1= 0.4595501 k=309, s1= 0.4595501: => new s1= 0.45951955 k=310, s1= 0.45951955: => new s1= 0.4595275 k=311, s1= 0.4595275: => new s1= 0.45960638 k=312, s1= 0.45960638: => new s1= 0.45959978 k=313, s1= 0.45959978: => new s1= 0.45960485 k=314, s1= 0.45960485: => new s1= 0.45960049 k=315, s1= 0.45960049: => new s1= 0.45963691 k=316, s1= 0.45963691: => new s1= 0.45962929 k=317, s1= 0.45962929: => new s1= 0.45955478 k=318, s1= 0.45955478: => new s1= 0.45965137 k=319, s1= 0.45965137: => new s1= 0.4597476 k=320, s1= 0.4597476: => new s1= 0.45975317 k=321, s1= 0.45975317: => new s1= 0.4597508 k=322, s1= 0.4597508: => new s1= 0.45975624 k=323, s1= 0.45975624: => new s1= 0.45975359 k=324, s1= 0.45975359: => new s1= 0.45975884 k=325, s1= 0.45975884: => new s1= 0.45975601 k=326, s1= 0.45975601: => new s1= 0.45976111 k=327, s1= 0.45976111: => new s1= 0.45975815 k=328, s1= 0.45975815: => new s1= 0.45976314 k=329, s1= 0.45976314: => new s1= 0.45976009 k=330, s1= 0.45976009: => new s1= 0.459765 k=331, s1= 0.459765: => new s1= 0.45976188 k=332, s1= 0.45976188: => new s1= 0.45976674 k=333, s1= 0.45976674: => new s1= 0.45976356 k=334, s1= 0.45976356: => new s1= 0.45976839 k=335, s1= 0.45976839: => new s1= 0.45976517 k=336, s1= 0.45976517: => new s1= 0.45976997 k=337, s1= 0.45976997: => new s1= 0.45976673 k=338, s1= 0.45976673: => new s1= 0.45977152 k=339, s1= 0.45977152: => new s1= 0.45976825 k=340, s1= 0.45976825: => new s1= 0.45977303 k=341, s1= 0.45977303: => new s1= 0.45976975 k=342, s1= 0.45976975: => new s1= 0.45977453 k=343, s1= 0.45977453: => new s1= 0.45977123 k=344, s1= 0.45977123: => new s1= 0.459776 k=345, s1= 0.459776: => new s1= 0.45977269 k=346, s1= 0.45977269: => new s1= 0.45977746 k=347, s1= 0.45977746: => new s1= 0.45977414 k=348, s1= 0.45977414: => new s1= 0.45977891 k=349, s1= 0.45977891: => new s1= 0.45977558 k=350, s1= 0.45977558: => new s1= 0.45978035 k=351, s1= 0.45978035: => new s1= 0.45977702 k=352, s1= 0.45977702: => new s1= 0.45978179 k=353, s1= 0.45978179: => new s1= 0.45977846 k=354, s1= 0.45977846: => new s1= 0.45978323 k=355, s1= 0.45978323: => new s1= 0.45977989 k=356, s1= 0.45977989: => new s1= 0.45978467 k=357, s1= 0.45978467: => new s1= 0.45978132 k=358, s1= 0.45978132: => new s1= 0.4597861 k=359, s1= 0.4597861: => new s1= 0.45978275 k=360, s1= 0.45978275: => new s1= 0.45978753 k=361, s1= 0.45978753: => new s1= 0.45978418 k=362, s1= 0.45978418: => new s1= 0.45978895 k=363, s1= 0.45978895: => new s1= 0.45978559 k=364, s1= 0.45978559: => new s1= 0.45979037 k=365, s1= 0.45979037: => new s1= 0.45978701 k=366, s1= 0.45978701: => new s1= 0.4597918 k=367, s1= 0.4597918: => new s1= 0.45978844 k=368, s1= 0.45978844: => new s1= 0.45979323 k=369, s1= 0.45979323: => new s1= 0.45978987 k=370, s1= 0.45978987: => new s1= 0.45979466 k=371, s1= 0.45979466: => new s1= 0.45979128 k=372, s1= 0.45979128: => new s1= 0.45979608 k=373, s1= 0.45979608: => new s1= 0.45979271 k=374, s1= 0.45979271: => new s1= 0.45979751 k=375, s1= 0.45979751: => new s1= 0.45979413 k=376, s1= 0.45979413: => new s1= 0.45979893 k=377, s1= 0.45979893: => new s1= 0.45979554 k=378, s1= 0.45979554: => new s1= 0.45980035 k=379, s1= 0.45980035: => new s1= 0.45979697 k=380, s1= 0.45979697: => new s1= 0.45980178 k=381, s1= 0.45980178: => new s1= 0.45979839 k=382, s1= 0.45979839: => new s1= 0.4598032 k=383, s1= 0.4598032: => new s1= 0.45979981 k=384, s1= 0.45979981: => new s1= 0.45980463 k=385, s1= 0.45980463: => new s1= 0.45980123 k=386, s1= 0.45980123: => new s1= 0.45980604 k=387, s1= 0.45980604: => new s1= 0.45980264 k=388, s1= 0.45980264: => new s1= 0.45980746 k=389, s1= 0.45980746: => new s1= 0.45980405 k=390, s1= 0.45980405: => new s1= 0.45980887 k=391, s1= 0.45980887: => new s1= 0.45980547 k=392, s1= 0.45980547: => new s1= 0.4598103 k=393, s1= 0.4598103: => new s1= 0.45980688 k=394, s1= 0.45980688: => new s1= 0.45981172 k=395, s1= 0.45981172: => new s1= 0.45980831 k=396, s1= 0.45980831: => new s1= 0.45981314 k=397, s1= 0.45981314: => new s1= 0.45980972 k=398, s1= 0.45980972: => new s1= 0.45981456 k=399, s1= 0.45981456: => new s1= 0.45981114 k=400, s1= 0.45981114: => new s1= 0.45981597 k=401, s1= 0.45981597: => new s1= 0.45981255 k=402, s1= 0.45981255: => new s1= 0.45981739 k=403, s1= 0.45981739: => new s1= 0.45981395 k=404, s1= 0.45981395: => new s1= 0.4598188 k=405, s1= 0.4598188: => new s1= 0.45981536 k=406, s1= 0.45981536: => new s1= 0.45982021 k=407, s1= 0.45982021: => new s1= 0.45981678 k=408, s1= 0.45981678: => new s1= 0.45982162 k=409, s1= 0.45982162: => new s1= 0.45981818 k=410, s1= 0.45981818: => new s1= 0.45982304 k=411, s1= 0.45982304: => new s1= 0.4598196 k=412, s1= 0.4598196: => new s1= 0.45982445 k=413, s1= 0.45982445: => new s1= 0.459821 k=414, s1= 0.459821: => new s1= 0.45982585 k=415, s1= 0.45982585: => new s1= 0.4598224 k=416, s1= 0.4598224: => new s1= 0.45982725 k=417, s1= 0.45982725: => new s1= 0.4598238 k=418, s1= 0.4598238: => new s1= 0.45982866 k=419, s1= 0.45982866: => new s1= 0.45982521 k=420, s1= 0.45982521: => new s1= 0.45983006 k=421, s1= 0.45983006: => new s1= 0.45982661 k=422, s1= 0.45982661: => new s1= 0.45983147 k=423, s1= 0.45983147: => new s1= 0.45982801 k=424, s1= 0.45982801: => new s1= 0.45983287 k=425, s1= 0.45983287: => new s1= 0.45982941 k=426, s1= 0.45982941: => new s1= 0.45983427 k=427, s1= 0.45983427: => new s1= 0.45983079 k=428, s1= 0.45983079: => new s1= 0.45983567 k=429, s1= 0.45983567: => new s1= 0.45983219 k=430, s1= 0.45983219: => new s1= 0.45983707 k=431, s1= 0.45983707: => new s1= 0.45983359 k=432, s1= 0.45983359: => new s1= 0.45983847 k=433, s1= 0.45983847: => new s1= 0.45983498 k=434, s1= 0.45983498: => new s1= 0.45983986 k=435, s1= 0.45983986: => new s1= 0.45983637 k=436, s1= 0.45983637: => new s1= 0.45984125 k=437, s1= 0.45984125: => new s1= 0.45983777 k=438, s1= 0.45983777: => new s1= 0.45984265 k=439, s1= 0.45984265: => new s1= 0.45983915 k=440, s1= 0.45983915: => new s1= 0.45984403 k=441, s1= 0.45984403: => new s1= 0.45984054 k=442, s1= 0.45984054: => new s1= 0.45984543 k=443, s1= 0.45984543: => new s1= 0.45984192 k=444, s1= 0.45984192: => new s1= 0.45984681 k=445, s1= 0.45984681: => new s1= 0.45984331 k=446, s1= 0.45984331: => new s1= 0.45984821 k=447, s1= 0.45984821: => new s1= 0.4598447 k=448, s1= 0.4598447: => new s1= 0.45984959 k=449, s1= 0.45984959: => new s1= 0.45984607 k=450, s1= 0.45984607: => new s1= 0.45985097 k=451, s1= 0.45985097: => new s1= 0.45984745 k=452, s1= 0.45984745: => new s1= 0.45985235 k=453, s1= 0.45985235: => new s1= 0.45984883 k=454, s1= 0.45984883: => new s1= 0.45985373 k=455, s1= 0.45985373: => new s1= 0.45985022 k=456, s1= 0.45985022: => new s1= 0.45985512 k=457, s1= 0.45985512: => new s1= 0.4598516 k=458, s1= 0.4598516: => new s1= 0.4598565 k=459, s1= 0.4598565: => new s1= 0.45985297 k=460, s1= 0.45985297: => new s1= 0.45985787 k=461, s1= 0.45985787: => new s1= 0.45985434 k=462, s1= 0.45985434: => new s1= 0.45985925 k=463, s1= 0.45985925: => new s1= 0.45985571 k=464, s1= 0.45985571: => new s1= 0.45986062 k=465, s1= 0.45986062: => new s1= 0.45985708 k=466, s1= 0.45985708: => new s1= 0.45986199 k=467, s1= 0.45986199: => new s1= 0.45985844 k=468, s1= 0.45985844: => new s1= 0.45986335 k=469, s1= 0.45986335: => new s1= 0.45985981 k=470, s1= 0.45985981: => new s1= 0.45986472 k=471, s1= 0.45986472: => new s1= 0.45986118 k=472, s1= 0.45986118: => new s1= 0.45986609 k=473, s1= 0.45986609: => new s1= 0.45986254 k=474, s1= 0.45986254: => new s1= 0.45986747 k=475, s1= 0.45986747: => new s1= 0.4598639 k=476, s1= 0.4598639: => new s1= 0.45986882 k=477, s1= 0.45986882: => new s1= 0.45986526 k=478, s1= 0.45986526: => new s1= 0.45987018 k=479, s1= 0.45987018: => new s1= 0.45986662 k=480, s1= 0.45986662: => new s1= 0.45987154 k=481, s1= 0.45987154: => new s1= 0.45986797 k=482, s1= 0.45986797: => new s1= 0.4598729 k=483, s1= 0.4598729: => new s1= 0.45986934 k=484, s1= 0.45986934: => new s1= 0.45987427 k=485, s1= 0.45987427: => new s1= 0.45987069 k=486, s1= 0.45987069: => new s1= 0.45987562 k=487, s1= 0.45987562: => new s1= 0.45987205 k=488, s1= 0.45987205: => new s1= 0.45987699 k=489, s1= 0.45987699: => new s1= 0.4598734 k=490, s1= 0.4598734: => new s1= 0.45987834 k=491, s1= 0.45987834: => new s1= 0.45987475 k=492, s1= 0.45987475: => new s1= 0.45987969 k=493, s1= 0.45987969: => new s1= 0.4598761 k=494, s1= 0.4598761: => new s1= 0.45988104 k=495, s1= 0.45988104: => new s1= 0.45987744 k=496, s1= 0.45987744: => new s1= 0.45988239 k=497, s1= 0.45988239: => new s1= 0.45987879 k=498, s1= 0.45987879: => new s1= 0.45988374 k=499, s1= 0.45988374: => new s1= 0.45988015 k=500, s1= 0.45988015: => new s1= 0.45988509 k=501, s1= 0.45988509: => new s1= 0.45988149 k=502, s1= 0.45988149: => new s1= 0.45988643 k=503, s1= 0.45988643: => new s1= 0.45988282 k=504, s1= 0.45988282: => new s1= 0.45988777 k=505, s1= 0.45988777: => new s1= 0.45988415 k=506, s1= 0.45988415: => new s1= 0.4598891 k=507, s1= 0.4598891: => new s1= 0.4598855 k=508, s1= 0.4598855: => new s1= 0.45989045 k=509, s1= 0.45989045: => new s1= 0.45988684 k=510, s1= 0.45988684: => new s1= 0.45989178 k=511, s1= 0.45989178: => new s1= 0.45988816 k=512, s1= 0.45988816: => new s1= 0.45989312 k=513, s1= 0.45989312: => new s1= 0.4598895 k=514, s1= 0.4598895: => new s1= 0.45989445 k=515, s1= 0.45989445: => new s1= 0.45989082 k=516, s1= 0.45989082: => new s1= 0.45989578 k=517, s1= 0.45989578: => new s1= 0.45989216 k=518, s1= 0.45989216: => new s1= 0.45989712 k=519, s1= 0.45989712: => new s1= 0.45989349 k=520, s1= 0.45989349: => new s1= 0.45989844 k=521, s1= 0.45989844: => new s1= 0.45989482 k=522, s1= 0.45989482: => new s1= 0.45989978 k=523, s1= 0.45989978: => new s1= 0.45989614 k=524, s1= 0.45989614: => new s1= 0.4599011 k=525, s1= 0.4599011: => new s1= 0.45989746 k=526, s1= 0.45989746: => new s1= 0.45990243 k=527, s1= 0.45990243: => new s1= 0.45989879 k=528, s1= 0.45989879: => new s1= 0.45990375 k=529, s1= 0.45990375: => new s1= 0.4599001 k=530, s1= 0.4599001: => new s1= 0.45990506 k=531, s1= 0.45990506: => new s1= 0.45990141 k=532, s1= 0.45990141: => new s1= 0.45990638 k=533, s1= 0.45990638: => new s1= 0.45990272 k=534, s1= 0.45990272: => new s1= 0.45990769 k=535, s1= 0.45990769: => new s1= 0.45990403 k=536, s1= 0.45990403: => new s1= 0.45990901 k=537, s1= 0.45990901: => new s1= 0.45990534 k=538, s1= 0.45990534: => new s1= 0.45991032 k=539, s1= 0.45991032: => new s1= 0.45990666 k=540, s1= 0.45990666: => new s1= 0.45991163 k=541, s1= 0.45991163: => new s1= 0.45990797 k=542, s1= 0.45990797: => new s1= 0.45991294 k=543, s1= 0.45991294: => new s1= 0.45990927 k=544, s1= 0.45990927: => new s1= 0.45991425 k=545, s1= 0.45991425: => new s1= 0.45991057 k=546, s1= 0.45991057: => new s1= 0.45991555 k=547, s1= 0.45991555: => new s1= 0.45991188 k=548, s1= 0.45991188: => new s1= 0.45991686 k=549, s1= 0.45991686: => new s1= 0.45991318 k=550, s1= 0.45991318: => new s1= 0.45991817 k=551, s1= 0.45991817: => new s1= 0.45991448 k=552, s1= 0.45991448: => new s1= 0.45991946 k=553, s1= 0.45991946: => new s1= 0.45991577 k=554, s1= 0.45991577: => new s1= 0.45992077 k=555, s1= 0.45992077: => new s1= 0.45991708 k=556, s1= 0.45991708: => new s1= 0.45992207 k=557, s1= 0.45992207: => new s1= 0.45991838 k=558, s1= 0.45991838: => new s1= 0.45992337 k=559, s1= 0.45992337: => new s1= 0.45991968 k=560, s1= 0.45991968: => new s1= 0.45992467 k=561, s1= 0.45992467: => new s1= 0.45992098 k=562, s1= 0.45992098: => new s1= 0.45992597 k=563, s1= 0.45992597: => new s1= 0.45992226 k=564, s1= 0.45992226: => new s1= 0.45992726 k=565, s1= 0.45992726: => new s1= 0.45992355 k=566, s1= 0.45992355: => new s1= 0.45992854 k=567, s1= 0.45992854: => new s1= 0.45992483 k=568, s1= 0.45992483: => new s1= 0.45992983 k=569, s1= 0.45992983: => new s1= 0.45992612 k=570, s1= 0.45992612: => new s1= 0.45993112 k=571, s1= 0.45993112: => new s1= 0.4599274 k=572, s1= 0.4599274: => new s1= 0.4599324 k=573, s1= 0.4599324: => new s1= 0.45992868 k=574, s1= 0.45992868: => new s1= 0.45993368 k=575, s1= 0.45993368: => new s1= 0.45992995 k=576, s1= 0.45992995: => new s1= 0.45993495 k=577, s1= 0.45993495: => new s1= 0.45993123 k=578, s1= 0.45993123: => new s1= 0.45993623 k=579, s1= 0.45993623: => new s1= 0.4599325 k=580, s1= 0.4599325: => new s1= 0.45993751 k=581, s1= 0.45993751: => new s1= 0.45993378 k=582, s1= 0.45993378: => new s1= 0.45993878 k=583, s1= 0.45993878: => new s1= 0.45993505 k=584, s1= 0.45993505: => new s1= 0.45994006 k=585, s1= 0.45994006: => new s1= 0.45993632 k=586, s1= 0.45993632: => new s1= 0.45994132 k=587, s1= 0.45994132: => new s1= 0.45993759 k=588, s1= 0.45993759: => new s1= 0.4599426 k=589, s1= 0.4599426: => new s1= 0.45993886 k=590, s1= 0.45993886: => new s1= 0.45994386 k=591, s1= 0.45994386: => new s1= 0.45994012 k=592, s1= 0.45994012: => new s1= 0.45994513 k=593, s1= 0.45994513: => new s1= 0.45994138 k=594, s1= 0.45994138: => new s1= 0.4599464 k=595, s1= 0.4599464: => new s1= 0.45994265 k=596, s1= 0.45994265: => new s1= 0.45994766 k=597, s1= 0.45994766: => new s1= 0.45994391 k=598, s1= 0.45994391: => new s1= 0.45994892 k=599, s1= 0.45994892: => new s1= 0.45994517 k=600, s1= 0.45994517: => new s1= 0.45995019 k=601, s1= 0.45995019: => new s1= 0.45994643 k=602, s1= 0.45994643: => new s1= 0.45995145 k=603, s1= 0.45995145: => new s1= 0.45994768 k=604, s1= 0.45994768: => new s1= 0.4599527 k=605, s1= 0.4599527: => new s1= 0.45994893 k=606, s1= 0.45994893: => new s1= 0.45995395 k=607, s1= 0.45995395: => new s1= 0.45995018 k=608, s1= 0.45995018: => new s1= 0.4599552 k=609, s1= 0.4599552: => new s1= 0.45995143 k=610, s1= 0.45995143: => new s1= 0.45995646 k=611, s1= 0.45995646: => new s1= 0.45995268 k=612, s1= 0.45995268: => new s1= 0.45995771 k=613, s1= 0.45995771: => new s1= 0.45995394 k=614, s1= 0.45995394: => new s1= 0.45995897 k=615, s1= 0.45995897: => new s1= 0.45995518 k=616, s1= 0.45995518: => new s1= 0.45996021 k=617, s1= 0.45996021: => new s1= 0.45995642 k=618, s1= 0.45995642: => new s1= 0.45996145 k=619, s1= 0.45996145: => new s1= 0.45995767 k=620, s1= 0.45995767: => new s1= 0.45996269 k=621, s1= 0.45996269: => new s1= 0.45995891 k=622, s1= 0.45995891: => new s1= 0.45996394 k=623, s1= 0.45996394: => new s1= 0.45996014 k=624, s1= 0.45996014: => new s1= 0.45996516 k=625, s1= 0.45996516: => new s1= 0.45996137 k=626, s1= 0.45996137: => new s1= 0.45996641 k=627, s1= 0.45996641: => new s1= 0.45996261 k=628, s1= 0.45996261: => new s1= 0.45996764 k=629, s1= 0.45996764: => new s1= 0.45996383 k=630, s1= 0.45996383: => new s1= 0.45996887 k=631, s1= 0.45996887: => new s1= 0.45996507 k=632, s1= 0.45996507: => new s1= 0.4599701 k=633, s1= 0.4599701: => new s1= 0.4599663 k=634, s1= 0.4599663: => new s1= 0.45997133 k=635, s1= 0.45997133: => new s1= 0.45996752 k=636, s1= 0.45996752: => new s1= 0.45997256 k=637, s1= 0.45997256: => new s1= 0.45996875 k=638, s1= 0.45996875: => new s1= 0.45997379 k=639, s1= 0.45997379: => new s1= 0.45996998 k=640, s1= 0.45996998: => new s1= 0.45997502 k=641, s1= 0.45997502: => new s1= 0.4599712 k=642, s1= 0.4599712: => new s1= 0.45997624 k=643, s1= 0.45997624: => new s1= 0.45997241 k=644, s1= 0.45997241: => new s1= 0.45997747 k=645, s1= 0.45997747: => new s1= 0.45997364 k=646, s1= 0.45997364: => new s1= 0.45997869 k=647, s1= 0.45997869: => new s1= 0.45997487 k=648, s1= 0.45997487: => new s1= 0.45997992 k=649, s1= 0.45997992: => new s1= 0.45997608 k=650, s1= 0.45997608: => new s1= 0.45998112 k=651, s1= 0.45998112: => new s1= 0.45997729 k=652, s1= 0.45997729: => new s1= 0.45998233 k=653, s1= 0.45998233: => new s1= 0.45997849 k=654, s1= 0.45997849: => new s1= 0.45998355 k=655, s1= 0.45998355: => new s1= 0.45997971 k=656, s1= 0.45997971: => new s1= 0.45998476 k=657, s1= 0.45998476: => new s1= 0.45998092 k=658, s1= 0.45998092: => new s1= 0.45998597 k=659, s1= 0.45998597: => new s1= 0.45998212 k=660, s1= 0.45998212: => new s1= 0.45998717 k=661, s1= 0.45998717: => new s1= 0.45998332 k=662, s1= 0.45998332: => new s1= 0.45998838 k=663, s1= 0.45998838: => new s1= 0.45998453 k=664, s1= 0.45998453: => new s1= 0.45998958 k=665, s1= 0.45998958: => new s1= 0.45998573 k=666, s1= 0.45998573: => new s1= 0.45999078 k=667, s1= 0.45999078: => new s1= 0.45998693 k=668, s1= 0.45998693: => new s1= 0.45999199 k=669, s1= 0.45999199: => new s1= 0.45998813 k=670, s1= 0.45998813: => new s1= 0.45999319 k=671, s1= 0.45999319: => new s1= 0.45998933 k=672, s1= 0.45998933: => new s1= 0.45999439 k=673, s1= 0.45999439: => new s1= 0.45999052 k=674, s1= 0.45999052: => new s1= 0.45999559 k=675, s1= 0.45999559: => new s1= 0.45999172 k=676, s1= 0.45999172: => new s1= 0.45999678 k=677, s1= 0.45999678: => new s1= 0.45999291 k=678, s1= 0.45999291: => new s1= 0.45999797 k=679, s1= 0.45999797: => new s1= 0.4599941 k=680, s1= 0.4599941: => new s1= 0.45999916 k=681, s1= 0.45999916: => new s1= 0.45999529 k=682, s1= 0.45999529: => new s1= 0.46000035 k=683, s1= 0.46000035: => new s1= 0.45999646 k=684, s1= 0.45999646: => new s1= 0.46000153 k=685, s1= 0.46000153: => new s1= 0.45999765 k=686, s1= 0.45999765: => new s1= 0.46000272 k=687, s1= 0.46000272: => new s1= 0.45999883 k=688, s1= 0.45999883: => new s1= 0.46000391 k=689, s1= 0.46000391: => new s1= 0.46000002 k=690, s1= 0.46000002: => new s1= 0.46000508 k=691, s1= 0.46000508: => new s1= 0.4600012 k=692, s1= 0.4600012: => new s1= 0.46000627 k=693, s1= 0.46000627: => new s1= 0.46000237 k=694, s1= 0.46000237: => new s1= 0.46000745 k=695, s1= 0.46000745: => new s1= 0.46000355 k=696, s1= 0.46000355: => new s1= 0.46000862 k=697, s1= 0.46000862: => new s1= 0.46000473 k=698, s1= 0.46000473: => new s1= 0.4600098 k=699, s1= 0.4600098: => new s1= 0.4600059 k=700, s1= 0.4600059: => new s1= 0.46001098 k=701, s1= 0.46001098: => new s1= 0.46000707 k=702, s1= 0.46000707: => new s1= 0.46001215 k=703, s1= 0.46001215: => new s1= 0.46000824 k=704, s1= 0.46000824: => new s1= 0.46001332 k=705, s1= 0.46001332: => new s1= 0.46000941 k=706, s1= 0.46000941: => new s1= 0.46001449 k=707, s1= 0.46001449: => new s1= 0.46001058 k=708, s1= 0.46001058: => new s1= 0.46001565 k=709, s1= 0.46001565: => new s1= 0.46001175 k=710, s1= 0.46001175: => new s1= 0.46001682 k=711, s1= 0.46001682: => new s1= 0.4600129 k=712, s1= 0.4600129: => new s1= 0.46001798 k=713, s1= 0.46001798: => new s1= 0.46001406 k=714, s1= 0.46001406: => new s1= 0.46001914 k=715, s1= 0.46001914: => new s1= 0.46001522 k=716, s1= 0.46001522: => new s1= 0.46002031 k=717, s1= 0.46002031: => new s1= 0.46001638 k=718, s1= 0.46001638: => new s1= 0.46002147 k=719, s1= 0.46002147: => new s1= 0.46001754 k=720, s1= 0.46001754: => new s1= 0.46002262 k=721, s1= 0.46002262: => new s1= 0.4600187 k=722, s1= 0.4600187: => new s1= 0.46002379 k=723, s1= 0.46002379: => new s1= 0.46001985 k=724, s1= 0.46001985: => new s1= 0.46002494 k=725, s1= 0.46002494: => new s1= 0.46002101 k=726, s1= 0.46002101: => new s1= 0.4600261 k=727, s1= 0.4600261: => new s1= 0.46002215 k=728, s1= 0.46002215: => new s1= 0.46002724 k=729, s1= 0.46002724: => new s1= 0.4600233 k=730, s1= 0.4600233: => new s1= 0.46002839 k=731, s1= 0.46002839: => new s1= 0.46002445 k=732, s1= 0.46002445: => new s1= 0.46002954 k=733, s1= 0.46002954: => new s1= 0.4600256 k=734, s1= 0.4600256: => new s1= 0.46003068 k=735, s1= 0.46003068: => new s1= 0.46002674 k=736, s1= 0.46002674: => new s1= 0.46003183 k=737, s1= 0.46003183: => new s1= 0.46002788 k=738, s1= 0.46002788: => new s1= 0.46003297 k=739, s1= 0.46003297: => new s1= 0.46002902 k=740, s1= 0.46002902: => new s1= 0.46003411 k=741, s1= 0.46003411: => new s1= 0.46003016 k=742, s1= 0.46003016: => new s1= 0.46003525 k=743, s1= 0.46003525: => new s1= 0.4600313 k=744, s1= 0.4600313: => new s1= 0.46003639 k=745, s1= 0.46003639: => new s1= 0.46003244 k=746, s1= 0.46003244: => new s1= 0.46003753 k=747, s1= 0.46003753: => new s1= 0.46003358 k=748, s1= 0.46003358: => new s1= 0.46003867 k=749, s1= 0.46003867: => new s1= 0.46003471 k=750, s1= 0.46003471: => new s1= 0.4600398 k=751, s1= 0.4600398: => new s1= 0.46003583 k=752, s1= 0.46003583: => new s1= 0.46004093 k=753, s1= 0.46004093: => new s1= 0.46003696 k=754, s1= 0.46003696: => new s1= 0.46004206 k=755, s1= 0.46004206: => new s1= 0.46003809 k=756, s1= 0.46003809: => new s1= 0.46004319 k=757, s1= 0.46004319: => new s1= 0.46003922 k=758, s1= 0.46003922: => new s1= 0.46004432 k=759, s1= 0.46004432: => new s1= 0.46004034 k=760, s1= 0.46004034: => new s1= 0.46004544 k=761, s1= 0.46004544: => new s1= 0.46004147 k=762, s1= 0.46004147: => new s1= 0.46004657 k=763, s1= 0.46004657: => new s1= 0.46004259 k=764, s1= 0.46004259: => new s1= 0.46004769 k=765, s1= 0.46004769: => new s1= 0.46004371 k=766, s1= 0.46004371: => new s1= 0.46004881 k=767, s1= 0.46004881: => new s1= 0.46004482 k=768, s1= 0.46004482: => new s1= 0.46004992 k=769, s1= 0.46004992: => new s1= 0.46004594 k=770, s1= 0.46004594: => new s1= 0.46005104 k=771, s1= 0.46005104: => new s1= 0.46004705 k=772, s1= 0.46004705: => new s1= 0.46005216 k=773, s1= 0.46005216: => new s1= 0.46004816 k=774, s1= 0.46004816: => new s1= 0.46005328 k=775, s1= 0.46005328: => new s1= 0.46004928 k=776, s1= 0.46004928: => new s1= 0.46005439 k=777, s1= 0.46005439: => new s1= 0.46005039 k=778, s1= 0.46005039: => new s1= 0.4600555 k=779, s1= 0.4600555: => new s1= 0.4600515 k=780, s1= 0.4600515: => new s1= 0.46005661 k=781, s1= 0.46005661: => new s1= 0.46005261 k=782, s1= 0.46005261: => new s1= 0.46005772 k=783, s1= 0.46005772: => new s1= 0.46005371 k=784, s1= 0.46005371: => new s1= 0.46005882 k=785, s1= 0.46005882: => new s1= 0.46005481 k=786, s1= 0.46005481: => new s1= 0.46005993 k=787, s1= 0.46005993: => new s1= 0.46005591 k=788, s1= 0.46005591: => new s1= 0.46006103 k=789, s1= 0.46006103: => new s1= 0.46005702 k=790, s1= 0.46005702: => new s1= 0.46006213 k=791, s1= 0.46006213: => new s1= 0.46005812 k=792, s1= 0.46005812: => new s1= 0.46006323 k=793, s1= 0.46006323: => new s1= 0.46005921 k=794, s1= 0.46005921: => new s1= 0.46006433 k=795, s1= 0.46006433: => new s1= 0.46006032 k=796, s1= 0.46006032: => new s1= 0.46006543 k=797, s1= 0.46006543: => new s1= 0.46006141 k=798, s1= 0.46006141: => new s1= 0.46006652 k=799, s1= 0.46006652: => new s1= 0.4600625 k=800, s1= 0.4600625: => new s1= 0.46006761 k=801, s1= 0.46006761: => new s1= 0.46006358 k=802, s1= 0.46006358: => new s1= 0.4600687 k=803, s1= 0.4600687: => new s1= 0.46006467 k=804, s1= 0.46006467: => new s1= 0.46006979 k=805, s1= 0.46006979: => new s1= 0.46006576 k=806, s1= 0.46006576: => new s1= 0.46007089 k=807, s1= 0.46007089: => new s1= 0.46006685 k=808, s1= 0.46006685: => new s1= 0.46007197 k=809, s1= 0.46007197: => new s1= 0.46006793 k=810, s1= 0.46006793: => new s1= 0.46007306 k=811, s1= 0.46007306: => new s1= 0.46006902 k=812, s1= 0.46006902: => new s1= 0.46007414 k=813, s1= 0.46007414: => new s1= 0.46007011 k=814, s1= 0.46007011: => new s1= 0.46007522 k=815, s1= 0.46007522: => new s1= 0.46007119 k=816, s1= 0.46007119: => new s1= 0.4600763 k=817, s1= 0.4600763: => new s1= 0.46007226 k=818, s1= 0.46007226: => new s1= 0.46007739 k=819, s1= 0.46007739: => new s1= 0.46007334 k=820, s1= 0.46007334: => new s1= 0.46007847 k=821, s1= 0.46007847: => new s1= 0.46007441 k=822, s1= 0.46007441: => new s1= 0.46007954 k=823, s1= 0.46007954: => new s1= 0.46007549 k=824, s1= 0.46007549: => new s1= 0.46008062 k=825, s1= 0.46008062: => new s1= 0.46007657 k=826, s1= 0.46007657: => new s1= 0.4600817 k=827, s1= 0.4600817: => new s1= 0.46007764 k=828, s1= 0.46007764: => new s1= 0.46008276 k=829, s1= 0.46008276: => new s1= 0.4600787 k=830, s1= 0.4600787: => new s1= 0.46008384 k=831, s1= 0.46008384: => new s1= 0.46007977 k=832, s1= 0.46007977: => new s1= 0.4600849 k=833, s1= 0.4600849: => new s1= 0.46008084 k=834, s1= 0.46008084: => new s1= 0.46008598 k=835, s1= 0.46008598: => new s1= 0.46008191 k=836, s1= 0.46008191: => new s1= 0.46008704 k=837, s1= 0.46008704: => new s1= 0.46008297 k=838, s1= 0.46008297: => new s1= 0.4600881 k=839, s1= 0.4600881: => new s1= 0.46008402 k=840, s1= 0.46008402: => new s1= 0.46008916 k=841, s1= 0.46008916: => new s1= 0.46008508 k=842, s1= 0.46008508: => new s1= 0.46009022 k=843, s1= 0.46009022: => new s1= 0.46008614 k=844, s1= 0.46008614: => new s1= 0.46009127 k=845, s1= 0.46009127: => new s1= 0.4600872 k=846, s1= 0.4600872: => new s1= 0.46009233 k=847, s1= 0.46009233: => new s1= 0.46008825 k=848, s1= 0.46008825: => new s1= 0.46009338 k=849, s1= 0.46009338: => new s1= 0.4600893 k=850, s1= 0.4600893: => new s1= 0.46009443 k=851, s1= 0.46009443: => new s1= 0.46009035 k=852, s1= 0.46009035: => new s1= 0.46009549 k=853, s1= 0.46009549: => new s1= 0.46009141 k=854, s1= 0.46009141: => new s1= 0.46009655 k=855, s1= 0.46009655: => new s1= 0.46009246 k=856, s1= 0.46009246: => new s1= 0.4600976 k=857, s1= 0.4600976: => new s1= 0.46009351 k=858, s1= 0.46009351: => new s1= 0.46009865 k=859, s1= 0.46009865: => new s1= 0.46009456 k=860, s1= 0.46009456: => new s1= 0.46009969 k=861, s1= 0.46009969: => new s1= 0.4600956 k=862, s1= 0.4600956: => new s1= 0.46010074 k=863, s1= 0.46010074: => new s1= 0.46009665 k=864, s1= 0.46009665: => new s1= 0.46010179 k=865, s1= 0.46010179: => new s1= 0.4600977 k=866, s1= 0.4600977: => new s1= 0.46010284 k=867, s1= 0.46010284: => new s1= 0.46009874 k=868, s1= 0.46009874: => new s1= 0.46010388 k=869, s1= 0.46010388: => new s1= 0.46009977 k=870, s1= 0.46009977: => new s1= 0.46010492 k=871, s1= 0.46010492: => new s1= 0.46010081 k=872, s1= 0.46010081: => new s1= 0.46010596 k=873, s1= 0.46010596: => new s1= 0.46010185 k=874, s1= 0.46010185: => new s1= 0.46010699 k=875, s1= 0.46010699: => new s1= 0.46010288 k=876, s1= 0.46010288: => new s1= 0.46010802 k=877, s1= 0.46010802: => new s1= 0.46010391 k=878, s1= 0.46010391: => new s1= 0.46010905 k=879, s1= 0.46010905: => new s1= 0.46010494 k=880, s1= 0.46010494: => new s1= 0.46011009 k=881, s1= 0.46011009: => new s1= 0.46010597 k=882, s1= 0.46010597: => new s1= 0.46011111 k=883, s1= 0.46011111: => new s1= 0.460107 k=884, s1= 0.460107: => new s1= 0.46011215 k=885, s1= 0.46011215: => new s1= 0.46010802 k=886, s1= 0.46010802: => new s1= 0.46011317 k=887, s1= 0.46011317: => new s1= 0.46010905 k=888, s1= 0.46010905: => new s1= 0.4601142 k=889, s1= 0.4601142: => new s1= 0.46011007 k=890, s1= 0.46011007: => new s1= 0.46011522 k=891, s1= 0.46011522: => new s1= 0.46011109 k=892, s1= 0.46011109: => new s1= 0.46011624 k=893, s1= 0.46011624: => new s1= 0.46011211 k=894, s1= 0.46011211: => new s1= 0.46011726 k=895, s1= 0.46011726: => new s1= 0.46011313 k=896, s1= 0.46011313: => new s1= 0.46011829 k=897, s1= 0.46011829: => new s1= 0.46011415 k=898, s1= 0.46011415: => new s1= 0.4601193 k=899, s1= 0.4601193: => new s1= 0.46011517 k=900, s1= 0.46011517: => new s1= 0.46012032 k=901, s1= 0.46012032: => new s1= 0.46011618 k=902, s1= 0.46011618: => new s1= 0.46012134 k=903, s1= 0.46012134: => new s1= 0.4601172 k=904, s1= 0.4601172: => new s1= 0.46012236 k=905, s1= 0.46012236: => new s1= 0.46011821 k=906, s1= 0.46011821: => new s1= 0.46012336 k=907, s1= 0.46012336: => new s1= 0.46011922 k=908, s1= 0.46011922: => new s1= 0.46012438 k=909, s1= 0.46012438: => new s1= 0.46012023 k=910, s1= 0.46012023: => new s1= 0.46012539 k=911, s1= 0.46012539: => new s1= 0.46012125 k=912, s1= 0.46012125: => new s1= 0.4601264 k=913, s1= 0.4601264: => new s1= 0.46012225 k=914, s1= 0.46012225: => new s1= 0.46012741 k=915, s1= 0.46012741: => new s1= 0.46012326 k=916, s1= 0.46012326: => new s1= 0.46012841 k=917, s1= 0.46012841: => new s1= 0.46012426 k=918, s1= 0.46012426: => new s1= 0.46012942 k=919, s1= 0.46012942: => new s1= 0.46012526 k=920, s1= 0.46012526: => new s1= 0.46013043 k=921, s1= 0.46013043: => new s1= 0.46012627 k=922, s1= 0.46012627: => new s1= 0.46013143 k=923, s1= 0.46013143: => new s1= 0.46012727 k=924, s1= 0.46012727: => new s1= 0.46013243 k=925, s1= 0.46013243: => new s1= 0.46012827 k=926, s1= 0.46012827: => new s1= 0.46013343 k=927, s1= 0.46013343: => new s1= 0.46012926 k=928, s1= 0.46012926: => new s1= 0.46013443 k=929, s1= 0.46013443: => new s1= 0.46013026 k=930, s1= 0.46013026: => new s1= 0.46013543 k=931, s1= 0.46013543: => new s1= 0.46013125 k=932, s1= 0.46013125: => new s1= 0.46013642 k=933, s1= 0.46013642: => new s1= 0.46013225 k=934, s1= 0.46013225: => new s1= 0.46013741 k=935, s1= 0.46013741: => new s1= 0.46013324 k=936, s1= 0.46013324: => new s1= 0.46013841 k=937, s1= 0.46013841: => new s1= 0.46013423 k=938, s1= 0.46013423: => new s1= 0.4601394 k=939, s1= 0.4601394: => new s1= 0.46013522 k=940, s1= 0.46013522: => new s1= 0.46014038 k=941, s1= 0.46014038: => new s1= 0.4601362 k=942, s1= 0.4601362: => new s1= 0.46014136 k=943, s1= 0.46014136: => new s1= 0.46013719 k=944, s1= 0.46013719: => new s1= 0.46014235 k=945, s1= 0.46014235: => new s1= 0.46013817 k=946, s1= 0.46013817: => new s1= 0.46014334 k=947, s1= 0.46014334: => new s1= 0.46013916 k=948, s1= 0.46013916: => new s1= 0.46014432 k=949, s1= 0.46014432: => new s1= 0.46014013 k=950, s1= 0.46014013: => new s1= 0.4601453 k=951, s1= 0.4601453: => new s1= 0.46014111 k=952, s1= 0.46014111: => new s1= 0.46014627 k=953, s1= 0.46014627: => new s1= 0.46014209 k=954, s1= 0.46014209: => new s1= 0.46014726 k=955, s1= 0.46014726: => new s1= 0.46014307 k=956, s1= 0.46014307: => new s1= 0.46014823 k=957, s1= 0.46014823: => new s1= 0.46014405 k=958, s1= 0.46014405: => new s1= 0.46014922 k=959, s1= 0.46014922: => new s1= 0.46014502 k=960, s1= 0.46014502: => new s1= 0.46015019 k=961, s1= 0.46015019: => new s1= 0.460146 k=962, s1= 0.460146: => new s1= 0.46015117 k=963, s1= 0.46015117: => new s1= 0.46014697 k=964, s1= 0.46014697: => new s1= 0.46015214 k=965, s1= 0.46015214: => new s1= 0.46014793 k=966, s1= 0.46014793: => new s1= 0.46015311 k=967, s1= 0.46015311: => new s1= 0.46014891 k=968, s1= 0.46014891: => new s1= 0.46015408 k=969, s1= 0.46015408: => new s1= 0.46014987 k=970, s1= 0.46014987: => new s1= 0.46015505 k=971, s1= 0.46015505: => new s1= 0.46015084 k=972, s1= 0.46015084: => new s1= 0.46015601 k=973, s1= 0.46015601: => new s1= 0.4601518 k=974, s1= 0.4601518: => new s1= 0.46015697 k=975, s1= 0.46015697: => new s1= 0.46015276 k=976, s1= 0.46015276: => new s1= 0.46015793 k=977, s1= 0.46015793: => new s1= 0.46015372 k=978, s1= 0.46015372: => new s1= 0.4601589 k=979, s1= 0.4601589: => new s1= 0.46015469 k=980, s1= 0.46015469: => new s1= 0.46015986 k=981, s1= 0.46015986: => new s1= 0.46015565 k=982, s1= 0.46015565: => new s1= 0.46016082 k=983, s1= 0.46016082: => new s1= 0.46015661 k=984, s1= 0.46015661: => new s1= 0.46016179 k=985, s1= 0.46016179: => new s1= 0.46015757 k=986, s1= 0.46015757: => new s1= 0.46016275 k=987, s1= 0.46016275: => new s1= 0.46015852 k=988, s1= 0.46015852: => new s1= 0.46016369 k=989, s1= 0.46016369: => new s1= 0.46016327 k=990, s1= 0.46016327: => new s1= 0.46016959 k=991, s1= 0.46016959: => new s1= 0.46016575 k=992, s1= 0.46016575: => new s1= 0.46016705 k=993, s1= 0.46016705: => new s1= 0.46016105 k=994, s1= 0.46016105: => new s1= 0.46016534 k=995, s1= 0.46016534: => new s1= 0.46016477 k=996, s1= 0.46016477: => new s1= 0.46016737 k=997, s1= 0.46016737: => new s1= 0.46016417 k=998, s1= 0.46016417: => new s1= 0.46016994 k=999, s1= 0.46016994: => new s1= 0.46016603 Warning message: In BYlogreg(x0 = X.fs, y = y.fs, initwml = FALSE, addIntercept = FALSE, : No convergence in 1000 steps. > signif( + rbind(ML = coef(m.fsML), QL =coef(m.fsQL), + WBY0=coef(m.fsWBY.), BY0=coef(m.fs.BY.), + WBY =coef(m.fsWBY ), BY =coef(m.fs.BY) + ) + , 4) (Intercept) tenancy suppl.income logInc ML 0.9264 -1.850 0.8961 -0.332800 QL 0.5894 -1.790 0.8167 -0.266900 WBY0 0.2243 -1.655 0.5899 -0.002904 BY0 -0.3407 -1.764 0.7767 -0.001492 WBY 0.2333 -1.721 0.6135 -0.002831 BY NA NA NA NA > > > if(FALSE) { + ## *scaling* of X ( ?? <==> ?? 'sigma1' ) ------------------ + + ## no "W" (Mahalanobis fail because of *singular* X): + m.fs.BY100 <- BYlogreg(x0=100*X.fs, initwml=FALSE, + y=y.fs, + addIntercept=FALSE, trace=TRUE, maxhalf=18) + ## ==> no convergence + + X1c <- cbind(1, 100*X.fs[,-1]) + m.fsWBY1c <- BYlogreg(x0=X1c, y=y.fs, + addIntercept=FALSE, trace=TRUE, maxhalf=18) + ## ==> illegal singularity$kind + + }## not yet > > ###-------- Gamma ------------ > > ## Realistic "data" {from help(glmrob)}: > mu <- c(122.131, 53.0979, 39.9039, 33.9232, 28.007, + 24.923, 21.5747, 19.6971, 18.4516) > ns.resid <- c(-0.0338228, 0.0923228, 0.0525284, 0.0317426, -0.035954, + 0.00308925, -0.026637, -0.0353932, -0.0244761) > Vmu <- c(14915.9, 2819.38, 1592.32, 1150.78, 784.39, + 621.156, 465.467, 387.978, 340.462) > Hp2 <- robustbase:::Huberprop2 > ## Hp2. <- robustbase:::Huberprop2. > > ## was: phis <- 2^(-70:-1) -- but that was *not* reliable (on 32-bit e.g.) > phis <- 2^(-42:-1) > H1 <- sapply(phis, function(phi) + Hp2(phi, ns.resid=ns.resid, mu=mu, Vmu=Vmu, tcc = 1.345)) > ## H2 <- sapply(phis, function(phi) > ## Hp2.(phi, ns.resid=ns.resid, mu=mu, Vmu=Vmu, tcc = 1.345)) > dput(signif(H1)) c(9.91741, 9.88674, 9.89438, 9.88674, 9.88961, 9.88961, 9.88961, 9.88984, 9.88973, 9.88964, 9.8897, 9.88975, 9.88976, 9.88975, 9.88974, 9.88974, 9.88974, 9.88974, 9.88974, 9.88974, 9.88974, 9.88974, 9.88975, 9.88975, 9.88975, 9.33161, 8.70618, 8.39347, 8.23714, 8.15902, 8.12006, 7.16275, 3.38703, -0.0879886, -2.3322, -4.16929, -5.26821, -5.80526, -6.04822, -6.11538, -6.02613, -5.66718 ) > H2 <- c(9.91741, + 9.88674, 9.89438, 9.88674, 9.88961, 9.88961, 9.88961, 9.88984, + 9.88973, 9.88964, 9.8897, 9.88975, 9.88976, 9.88975, 9.88974, + 9.88974, 9.88974, 9.88974, 9.88974, 9.88974, 9.88974, 9.88974, + 9.88975, 9.88975, 9.88975, 9.33161, 8.70618, 8.39347, 8.23714, + 8.15902, 8.12006, 7.16275, 3.38703, -0.0879886, -2.3322, -4.16929, + -5.26821, -5.80526, -6.04822, -6.11538, -6.02613, -5.66718) > all.equal(H1,H2, tolerance = 0) # -> see 8.869e-7 [1] "Mean relative difference: 3.774807e-07" > stopifnot(all.equal(H1,H2, tolerance = 1e-5)) > > if(dev.interactive(TRUE)) # shows that phi < 1e-12 is doubtful + matplot(phis, cbind(H1,H2), log="x", ylim = rrange(H1), type="o") > > > proc.time() user system elapsed 10.75 1.17 11.89