1### This is originally from the R package 2#### 3#### rrcov : Scalable Robust Estimators with High Breakdown Point 4#### 5#### by Valentin Todorov 6 7## I would like to thank Peter Rousseeuw and Katrien van Driessen for 8## providing the initial code of this function. 9 10## This program is free software; you can redistribute it and/or modify 11## it under the terms of the GNU General Public License as published by 12## the Free Software Foundation; either version 2 of the License, or 13## (at your option) any later version. 14## 15## This program is distributed in the hope that it will be useful, 16## but WITHOUT ANY WARRANTY; without even the implied warranty of 17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18## GNU General Public License for more details. 19## 20## You should have received a copy of the GNU General Public License 21## along with this program; if not, a copy is available at 22## http://www.r-project.org/Licenses/ 23 24## No longer hidden in namespace : 25## easier to explain when user-available & documented if 26h.alpha.n <- function(alpha, n, p) { 27 ## Compute h(alpha) := size of subsample, given alpha, (n,p) 28 ## Same function for covMcd() and ltsReg() 29 n2 <- (n+p+1) %/% 2 30 floor(2 * n2 - n + 2 * (n - n2) * alpha) 31} 32 33## MM: the way it's set up, *must* be kept in sync with rrcov.control()'s 34## defaults --> ./rrcov.control.R : 35covMcd <- function(x, 36 cor = FALSE, 37 raw.only = FALSE, 38 alpha = control$ alpha, 39 nsamp = control$ nsamp, 40 nmini = control$ nmini, kmini = control$ kmini, 41 scalefn=control$scalefn, maxcsteps=control$maxcsteps, 42 initHsets = NULL, save.hsets = FALSE, names = TRUE, 43 seed = control$ seed, 44 tolSolve = control$ tolSolve, # had 1e-10 hardwired {now 1e-14 default} 45 trace = control$ trace, 46 use.correction = control$ use.correction, 47 wgtFUN = control$ wgtFUN, 48 control = rrcov.control()) 49{ 50 logdet.Lrg <- 50 ## <-- FIXME add to rrcov.control() and then use that 51 ## Analyze and validate the input parameters ... 52 if(length(seed) > 0) { 53 if(length(seed) < 3 || seed[1L] < 100) 54 stop("invalid 'seed'. Must be compatible with .Random.seed !") 55 if(exists(".Random.seed", envir=.GlobalEnv, inherits=FALSE)) { 56 seed.keep <- get(".Random.seed", envir=.GlobalEnv, inherits=FALSE) 57 on.exit(assign(".Random.seed", seed.keep, envir=.GlobalEnv)) 58 } 59 assign(".Random.seed", seed, envir=.GlobalEnv) 60 } 61 62 ## For back compatibility, as some new args did not exist pre 2013-04, 63 ## and callers of covMcd() may use a "too small" 'control' list: 64 defCtrl <- if(missing(control)) control else rrcov.control() 65 if(missing(wgtFUN)) getDefCtrl("wgtFUN", defCtrl) 66 if(is.null (nmini)) getDefCtrl("nmini", defCtrl) 67 68 ## vt::03.02.2006 - added options "best" and "exact" for nsamp 69 ## nsamp will be further analized in the wrapper .fastmcd() 70 if(is.numeric(nsamp) && nsamp <= 0) 71 stop("Invalid number of trials nsamp = ",nsamp, "!") 72 73 if(is.data.frame(x)) 74 x <- data.matrix(x, rownames.force=FALSE) 75 else if (!is.matrix(x)) 76 x <- matrix(x, length(x), 1, 77 dimnames = if(names) list(names(x), deparse(substitute(x)))) 78 79 if(!names) dimnames(x) <- NULL # (speedup) 80 ## drop all rows with missing values (!!) : 81 ok <- is.finite(x %*% rep.int(1, ncol(x))) 82 x <- x[ok, , drop = FALSE] 83 if(!length(dx <- dim(x))) 84 stop("All observations have missing values!") 85 n <- dx[1]; p <- dx[2] 86 if(names) dimn <- dimnames(x) 87 ## h(alpha) , the size of the subsamples 88 h <- h.alpha.n(alpha, n, p) 89 if(n <= p + 1) # ==> floor((n+p+1)/2) > n - 1 -- not Ok 90 stop(if (n <= p) # absolute barrier! 91 "n <= p -- you can't be serious!" 92 else "n == p+1 is too small sample size for MCD") 93 ## else 94 if(n < 2 * p) { ## p+1 < n < 2p 95 warning("n < 2 * p, i.e., possibly too small sample size") 96 ## was stop("Need at least 2*(number of variables) observations ") 97 } 98 ## jmin <- (n + p + 1) %/% 2 99 ## if(alpha < 1/2) ## FIXME? shouldn't we rather test 'alpha < jmin/n' ? 100 ## stop("The MCD must cover at least", jmin, "observations") 101 ## MM: I think this should be sufficient; 102 ## we should even omit the (n < 2p) warning 103 if(h > n) 104 stop("Sample size n < h(alpha; n,p) := size of \"good\" subsample") 105 else if(2*h < n) 106 warning("subsample size h < n/2 may be too small") 107 108 if(is.character(wgtFUN)) { 109 if(is.function(mkWfun <- .wgtFUN.covMcd[[wgtFUN]])) 110 wgtFUN <- mkWfun(p=p, n=n, control) 111 } 112 if(!is.function(wgtFUN)) 113 stop(gettextf("'wgtFUN' must be a function or one of the strings %s.", 114 pasteK(paste0('"',names(.wgtFUN.covMcd),'"'))), domain=NA) 115 116 ## vt::03.02.2006 - raw.cnp2 and cnp2 are vectors of size 2 and will 117 ## contain the correction factors (concistency and finite sample) 118 ## for the raw and reweighted estimates respectively. Set them 119 ## initially to 1. If use.correction is false (not the default), 120 ## the finite sample correction factor will not be used 121 ## (neither for the raw estimates nor for the reweighted ones) 122 raw.cnp2 <- cnp2 <- c(1,1) 123 124 ans <- list(call = match.call(), nsamp = nsamp, 125 method = sprintf("MCD(alpha=%g ==> h=%d)", alpha, h)) 126 127 if(h == n) { ## <==> alpha ~= 1 : Just compute the classical estimates -------- 128 mcd <- cov(x) #MM: was cov.wt(x)$cov 129 loc <- as.vector(colMeans(x)) 130 obj <- determinant(mcd, logarithm = TRUE)$modulus[1] 131 if ( -obj/p > logdet.Lrg ) { 132 ans$cov <- mcd 133 if(names) dimnames(ans$cov) <- list(dimn[[2]], dimn[[2]]) 134 if (cor) 135 ans$cor <- cov2cor(ans$cov) 136 ans$center <- loc 137 if(names && length(dimn[[2]])) 138 names(ans$center) <- dimn[[2]] 139 ans$n.obs <- n 140 ans$singularity <- list(kind = "classical") 141 weights <- 1 142 } 143 else { 144 mah <- mahalanobis(x, loc, mcd, tol = tolSolve) 145 ## VT:: 01.09.2004 - bug in alpha=1 146 weights <- wgtFUN(mah) # 0/1 147 sum.w <- sum(weights) 148 ans <- c(ans, cov.wt(x, wt = weights, cor = cor)) 149 ## cov.wt() -> list("cov", "center", "n.obs", ["wt", "cor"]) 150 ## Consistency factor for reweighted MCD -- ok for default wgtFUN only: FIXME 151 if(sum.w != n) { 152 cnp2[1] <- .MCDcons(p, sum.w/n) 153 ans$cov <- ans$cov * cnp2[1] 154 } 155 obj <- determinant(mcd, logarithm = TRUE)$modulus[1] 156 if( -obj/p > logdet.Lrg ) { 157 ans$singularity <- list(kind = "reweighted.MCD") 158 } 159 else { 160 mah <- mahalanobis(x, ans$center, ans$cov, tol = tolSolve) 161 weights <- wgtFUN(mah) # 0/1 162 } 163 } 164 165 ans$alpha <- alpha 166 ans$quan <- h 167 ans$raw.cov <- mcd 168 ans$raw.center <- loc 169 if(names && !is.null(nms <- dimn[[2]])) { 170 names(ans$raw.center) <- nms 171 dimnames(ans$raw.cov) <- list(nms,nms) 172 } 173 ans$crit <- obj # was exp(obj); but log-scale is "robust" against under/overflow 174 ans$method <- paste(ans$method, 175 "\nalpha = 1: The minimum covariance determinant estimates based on", 176 n, "observations \nare equal to the classical estimates.") 177 ans$mcd.wt <- rep.int(NA, length(ok)) 178 ans$mcd.wt[ok] <- weights 179 if(names && length(dimn[[1]])) 180 names(ans$mcd.wt) <- dimn[[1]] 181 ans$wt <- NULL 182 ans$X <- x 183 if(names) { 184 if(length(dimn[[1]])) 185 dimnames(ans$X)[[1]] <- names(ans$mcd.wt)[ok] 186 else 187 dimnames(ans$X) <- list(seq(along = ok)[ok], NULL) 188 } 189 if(trace) 190 cat(ans$method, "\n") 191 ans$raw.cnp2 <- raw.cnp2 192 ans$cnp2 <- cnp2 193 class(ans) <- "mcd" 194 return(ans) 195 } ## end { alpha = 1 <==> h = n } 196 197 mcd <- if(nsamp == "deterministic") { 198 ans$method <- paste("Deterministic", ans$method) 199 .detmcd (x, h, hsets.init = initHsets, 200 save.hsets=save.hsets, # full.h=full.h, 201 scalefn=scalefn, maxcsteps=maxcsteps, trace=as.integer(trace), 202 names=names) 203 } else { 204 ans$method <- paste0("Fast ", ans$method, "; nsamp = ", nsamp, 205 "; (n,k)mini = (", nmini,",",kmini,")") 206 .fastmcd(x, h, nsamp, nmini, kmini, trace=as.integer(trace)) 207 } 208 209 ## Compute the consistency correction factor for the raw MCD 210 ## (see calfa in Croux and Haesbroeck) 211 calpha <- .MCDcons(p, h/n) ## VT::19.3.2007 212 correct <- if(use.correction) .MCDcnp2(p, n, alpha) else 1. 213 raw.cnp2 <- c(calpha, correct) 214 215 if(p == 1) { 216 ## ==> Compute univariate location and scale estimates 217 ans$method <- paste("Univariate", ans$method) 218 scale <- sqrt(calpha * correct) * as.double(mcd$initcovariance) 219 center <- as.double(mcd$initmean) 220 if(abs(scale - 0) < 1e-07) { 221 ans$singularity <- list(kind = "identicalObs", q = h) 222 ans$raw.cov <- ans$cov <- matrix(0) 223 ans$raw.center <- ans$center <- center 224 ans$n.obs <- n 225 ans$alpha <- alpha 226 ans$quan <- h 227 if(names && !is.null(nms <- dimn[[2]][1])) { 228 names(ans$raw.center) <- names(ans$center) <- nms 229 dimnames(ans$raw.cov) <- dimnames(ans$cov) <- list(nms,nms) 230 } 231 ans$crit <- -Inf # = log(0) 232 weights <- as.numeric(abs(x - center) < 1e-07) # 0 / 1 233 } ## end { scale ~= 0 } 234 else { 235 ## Compute the weights for the raw MCD in case p=1 236 weights <- wgtFUN(((x - center)/scale)^2) # 0/1 237 sum.w <- sum(weights) 238 ans <- c(ans, cov.wt(x, wt = weights, cor=cor)) 239 240 if(sum.w != n) { 241 cdelta.rew <- .MCDcons(p, sum.w/n) ## VT::19.3.2007 242 correct.rew <- if(use.correction) .MCDcnp2.rew(p, n, alpha) else 1. 243 cnp2 <- c(cdelta.rew, correct.rew) 244 ans$cov <- cdelta.rew * correct.rew * ans$cov 245 } 246 ans$alpha <- alpha 247 ans$quan <- h 248 ans$raw.cov <- as.matrix(scale^2) 249 ans$raw.center <- as.vector(center) 250 if(names && !is.null(nms <- dimn[[2]][1])) { 251 dimnames(ans$raw.cov) <- list(nms,nms) 252 names(ans$raw.center) <- nms 253 } 254 ans$crit <- ## log(det) = 255 log(sum(sort((x - as.double(mcd$initmean))^2, partial = h)[1:h])/max(1,h-1)) 256 center <- ans$center 257 scale <- as.vector(sqrt(ans$cov)) 258 weights <- wgtFUN(((x - center)/scale)^2) 259 } ## end{ scale > 0 } 260 } ## end p=1 261 262 else { ## p >= 2 : --------------------------------------------------------- 263 264 ## Apply correction factor to the raw estimates 265 ## and use them to compute weights 266 mcd$initcovariance <- matrix(calpha * correct * mcd$initcovariance, p,p) 267 if(raw.only || mcd$exactfit != 0) { 268 ## If not all observations are in general position, i.e. more than 269 ## h observations lie on a hyperplane, the program still yields 270 ## the MCD location and scatter matrix, the latter being singular 271 ## (as it should be), as well as the equation of the hyperplane. 272 273 dim(mcd$coeff) <- c(5, p) 274 ans$cov <- ans$raw.cov <- mcd$initcovariance 275 ans$center <- ans$raw.center <- as.vector(mcd$initmean) 276 277 if(names && !is.null(nms <- dimn[[2]])) { 278 dimnames(ans$cov) <- list(nms, nms) 279 names(ans$center) <- nms 280 } 281 ans$n.obs <- n 282 283 if(raw.only) { 284 ans$raw.only <- TRUE 285 } else { 286 ## no longer relevant: 287 ## if(mcd$exactfit == -1) 288 ## stop("The program allows for at most ", mcd$kount, " observations.") 289 ## if(mcd$exactfit == -2) 290 ## stop("The program allows for at most ", mcd$kount, " variables.") 291 if(!(mcd$exactfit %in% c(1,2,3))) 292 stop("Unexpected 'exactfit' code ", mcd$exactfit, ". Please report!") 293 ## new (2007-01) and *instead* of older long 'method' extension; 294 ## the old message is still *printed* via .MCDsingularityMsg() 295 ## 296 ## exactfit is now *passed* to result instead of coded into 'message': 297 ans$singularity <- 298 list(kind = "on.hyperplane", exactCode = mcd$exactfit, 299 p = p, h = h, count = mcd$kount, coeff = mcd$coeff[1,]) 300 } 301 ans$alpha <- alpha 302 ans$quan <- h 303 if(names && !is.null(nms <- dimn[[2]])) { 304 names(ans$raw.center) <- nms 305 dimnames(ans$raw.cov) <- list(nms,nms) 306 } 307 ans$crit <- -Inf # = log(0) 308 weights <- mcd$weights 309 310 } ## end (raw.only || exact fit) 311 312 else { ## have general position (exactfit == 0) : ------------------------ 313 314 ## FIXME? here, we assume that mcd$initcovariance is not singular: 315 mah <- mahalanobis(x, mcd$initmean, mcd$initcovariance, tol = tolSolve) 316 weights <- wgtFUN(mah) 317 sum.w <- sum(weights) 318 ans <- c(ans, cov.wt(x, wt = weights, cor=cor)) 319 ## simple check for singularity, much cheaper than determinant() below: 320 sing.rewt <- any(apply(ans$cov == 0, 2, all)) 321 322 ## Compute and apply the consistency correction factor for 323 ## the reweighted cov 324 if(!sing.rewt && sum.w != n) { 325 cdelta.rew <- .MCDcons(p, sum.w/n) ## VT::19.3.2007 326 correct.rew <- if(use.correction) .MCDcnp2.rew(p, n, alpha) else 1. 327 cnp2 <- c(cdelta.rew, correct.rew) 328 ans$cov <- cdelta.rew * correct.rew * ans$cov 329 } 330 331 ##vt:: add also the best found subsample to the result list 332 ans$best <- sort(as.vector(mcd$best)) 333 334 ans$alpha <- alpha 335 ans$quan <- h 336 ans$raw.cov <- mcd$initcovariance 337 ans$raw.center <- as.vector(mcd$initmean) 338 if(names && !is.null(nms <- dimn[[2]])) { 339 names(ans$raw.center) <- nms 340 dimnames(ans$raw.cov) <- list(nms,nms) 341 } 342 ans$raw.weights <- weights 343 ans$crit <- mcd$mcdestimate # now in log scale! 344 ## 'mah' already computed above 345 ans$raw.mah <- mah ## mahalanobis(x, ans$raw.center, ans$raw.cov, tol = tolSolve) 346 ## Check if the reweighted scatter matrix is singular. 347 if(sing.rewt || - determinant(ans$cov, logarithm = TRUE)$modulus[1]/p > logdet.Lrg) { 348 ans$singularity <- list(kind = paste0("reweighted.MCD", 349 if(sing.rewt)"(zero col.)")) 350 ans$mah <- mah 351 } 352 else { 353 mah <- mahalanobis(x, ans$center, ans$cov, tol = tolSolve) 354 ans$mah <- mah 355 weights <- wgtFUN(mah) 356 } 357 } ## end{ not exact fit } 358 359 } ## end{ p >= 2 } 360 361 ans$mcd.wt <- rep.int(NA, length(ok)) 362 ans$mcd.wt[ok] <- weights 363 if(names) { 364 if(length(dimn[[1]])) 365 names(ans$mcd.wt) <- dimn[[1]] 366 if(length(dimn[[1]])) 367 dimnames(x)[[1]] <- names(ans$mcd.wt)[ok] 368 else 369 dimnames(x) <- list(seq(along = ok)[ok], NULL) 370 } 371 ans$X <- x 372 ans$wt <- NULL 373 if(trace) 374 cat(ans$method, "\n") 375 ans$raw.cnp2 <- raw.cnp2 376 ans$cnp2 <- cnp2 377 if(nsamp == "deterministic") 378 ans <- c(ans, mcd[c("iBest","n.csteps", if(save.hsets) "initHsets")]) 379 class(ans) <- "mcd" 380 ## warn if we have a singularity: 381 if(is.list(ans$singularity)) 382 warning(paste(strwrap(.MCDsingularityMsg(ans$singularity, ans$n.obs)), collapse="\n"), 383 domain=NA) 384 ## return 385 ans 386} ## {covMcd} 387 388smoothWgt <- function(x, c, h) { 389 ## currently drops all attributes, dim(), names(), etc 390 ## maybe add 'keep.attributes = FALSE' (and pass to C) 391 .Call(R_wgt_flex, x, c, h) 392} 393 394##' Martin Maechler's simple proposal for an *adaptive* cutoff 395##' i.e., one which does *not* reject outliers in good samples asymptotically 396.MCDadaptWgt.c <- function(n,p) { 397 eps <- 0.4 / n ^ 0.6 # => 1-eps(n=100) ~= 0.975; 1-eps(n=10) ~= 0.90 398 ## using upper tail: 399 qchisq(eps, p, lower.tail=FALSE) 400} 401 402 403## Default wgtFUN()s : 404.wgtFUN.covMcd <- 405 list("01.original" = function(p, ...) { 406 cMah <- qchisq(0.975, p) 407 function(d) as.numeric(d < cMah) 408 }, 409 "01.flex" = function(p, n, control) { ## 'control$beta' instead of 0.975 410 ## FIXME: update rrcov.control() to accept 'beta' 411 stopifnot(is.1num(beta <- control$beta), 0 <= beta, beta <= 1) 412 cMah <- qchisq(beta, p) 413 function(d) as.numeric(d < cMah) 414 }, 415 "01.adaptive" = function(p, n, ...) { ## 'beta_n' instead of 0.975 416 cMah <- .MCDadaptWgt.c(n,p) 417 function(d) as.numeric(d < cMah) 418 }, 419 "sm1.orig" = function(p, n, ...) { 420 cMah <- qchisq(0.975, p) 421 function(d) smoothWgt(d, c = cMah, h = 1) 422 }, 423 "sm2.orig" = function(p, n, ...) { 424 cMah <- qchisq(0.975, p) 425 function(d) smoothWgt(d, c = cMah, h = 2) 426 }, 427 "sm1.adaptive" = function(p, n, ...) { 428 cMah <- .MCDadaptWgt.c(n,p) 429 function(d) smoothWgt(d, c = cMah, h = 1) 430 }, 431 "sm2.adaptive" = function(p, n, ...) { 432 cMah <- .MCDadaptWgt.c(n,p) 433 function(d) smoothWgt(d, c = cMah, h = 2) 434 }, 435 "smE.adaptive" = function(p, n, ...) { 436 cMah <- .MCDadaptWgt.c(n,p) 437 ## TODO: find "theory" for h = f(cMah), or better c=f1(n,p); h=f2(n,p) 438 function(d) smoothWgt(d, c = cMah, h = max(2, cMah/4)) 439 } 440 ) 441 442 443.MCDsingularityMsg <- function(singList, n.obs) 444{ 445 stopifnot(is.list(singList)) 446 447 switch(singList$kind, 448 "classical" = { 449 "The classical covariance matrix is singular." 450 }, 451 "reweighted.MCD" = { 452 "The reweighted MCD scatter matrix is singular." 453 }, 454 "identicalObs" = { 455 sprintf("Initial scale 0 because more than 'h' (=%d) observations are identical.", 456 singList$q) 457 }, 458 "on.hyperplane" = { 459 stopifnot(c("p", "count", "coeff") %in% names(singList)) 460 obsMsg <- function(m, n) 461 paste("There are", m, 462 "observations (in the entire dataset of", n, "obs.)", 463 "lying on the") 464 invisible(obsMsg)# <- codetools 465 with(singList, 466 c(switch(exactCode, 467 ## exactfit == 1 : 468 "The covariance matrix of the data is singular.", 469 ## exactfit == 2 : 470 c("The covariance matrix has become singular during", 471 "the iterations of the MCD algorithm."), 472 ## exactfit == 3: 473 paste0("The ", h, 474 "-th order statistic of the absolute deviation of variable ", 475 which(coeff == 1), " is zero.")), 476 477 if(p == 2) { 478 paste(obsMsg(count, n.obs), "line with equation ", 479 signif(coeff[1], digits= 5), "(x_i1-m_1) +", 480 signif(coeff[2], digits= 5), "(x_i2-m_2) = 0", 481 "with (m_1,m_2) the mean of these observations.") 482 } 483 else if(p == 3) { 484 paste(obsMsg(count, n.obs), "plane with equation ", 485 signif(coeff[1], digits= 5), "(x_i1-m_1) +", 486 signif(coeff[2], digits= 5), "(x_i2-m_2) +", 487 signif(coeff[3], digits= 5), "(x_i3-m_3) = 0", 488 "with (m_1,m_2) the mean of these observations." 489 ) 490 } 491 else { ## p > 3 ----------- 492 paste(obsMsg(count, n.obs), "hyperplane with equation ", 493 "a_1*(x_i1 - m_1) + ... + a_p*(x_ip - m_p) = 0", 494 "with (m_1, ..., m_p) the mean of these observations", 495 "and coefficients a_i from the vector a <- ", 496 paste(deparse(zapsmall(coeff)), collapse="\n ")) 497 })) 498 }, 499 ## Otherwise 500 stop("illegal 'singularity$kind'") 501 ) ## end{switch} 502} 503 504nobs.mcd <- function (object, ...) object$n.obs 505 506print.mcd <- function(x, digits = max(3, getOption("digits") - 3), print.gap = 2, ...) 507{ 508 cat("Minimum Covariance Determinant (MCD) estimator approximation.\n", 509 "Method: ", x$method, "\n", sep="") 510 if(!is.null(cl <- x$call)) { 511 cat("Call:\n") 512 dput(cl) 513 } 514 if(is.list(x$singularity)) 515 cat(strwrap(.MCDsingularityMsg(x$singularity, x$n.obs)), sep ="\n") 516 517 if(identical(x$nsamp, "deterministic")) 518 cat("iBest: ", pasteK(x$iBest), "; C-step iterations: ", pasteK(x$n.csteps), 519 "\n", sep="") 520 ## VT::29.03.2007 - solve a conflict with fastmcd() in package robust - 521 ## also returning an object of class "mcd" 522 xx <- NA 523 if(!is.null(x$crit)) 524 xx <- format(x$crit, digits = digits) 525 else if (!is.null(x$raw.objective)) 526 xx <- format(log(x$raw.objective), digits = digits) 527 cat("Log(Det.): ", xx , "\n\nRobust Estimate of Location:\n") 528 print(x$center, digits = digits, print.gap = print.gap, ...) 529 cat("Robust Estimate of Covariance:\n") 530 print(x$cov, digits = digits, print.gap = print.gap, ...) 531 invisible(x) 532} 533 534summary.mcd <- function(object, ...) 535{ 536 class(object) <- c("summary.mcd", class(object)) 537 object 538} 539 540print.summary.mcd <- 541 function(x, digits = max(3, getOption("digits") - 3), print.gap = 2, ...) 542{ 543 print.mcd(x, digits = digits, print.gap = print.gap, ...) # see above 544 545 ## hmm, maybe not *such* a good idea : 546 if(!is.null(x$cor)) { 547 cat("\nRobust Estimate of Correlation: \n") 548 dimnames(x$cor) <- dimnames(x$cov) 549 print(x$cor, digits = digits, print.gap = print.gap, ...) 550 } 551 552 cat("\nEigenvalues:\n") 553 print(eigen(x$cov, only.values = TRUE)$values, digits = digits, ...) 554 555 if(!is.null(x$mah)) { 556 cat("\nRobust Distances: \n") 557 print(summary(x$mah, digits = digits), digits = digits, ...) 558 } 559 if(!is.null(wt <- x$mcd.wt)) 560 summarizeRobWeights(wt, digits = digits) 561 invisible(x) 562} 563 564## NOTE: plot.mcd() is in ./covPlot.R ! 565## ---- ~~~~~~~~~~~ 566 567### Consistency and Finite Sample Correction Factors 568### .MCDcons() .MCDcnp2() & .MCDcnp2.rew() 569 570### now exported and documented in ../man/covMcd.Rd 571 572##' Compute the consistency correction factor for the MCD estimate 573##' (see calfa in Croux and Haesbroeck) 574##' @param p 575##' @param alpha alpha ~= h/n = quan/n 576##' also use for the reweighted MCD, calling with alpha = 'sum(weights)/n' 577MCDcons <- # <- *not* exported, but currently used in pkgs rrcov, rrcovNA 578.MCDcons <- function(p, alpha) 579{ 580 qalpha <- qchisq(alpha, p) 581 caI <- pgamma(qalpha/2, p/2 + 1) / alpha 582 1/caI 583} 584 585MCDcnp2 <- # <- *not* exported, but currently used in pkg rrcovNA 586##' Finite sample correction factor for raw MCD: 587.MCDcnp2 <- function(p, n, alpha) 588{ 589 stopifnot(0 <= alpha, alpha <= 1, length(alpha) == 1) 590 591 if(p > 2) { 592 ## "alfaq" "betaq" "qwaarden" 593 coeffqpkwad875 <- matrix(c(-0.455179464070565, 1.11192541278794, 2, 594 -0.294241208320834, 1.09649329149811, 3), 595 ncol = 2) 596 coeffqpkwad500 <- matrix(c(-1.42764571687802, 1.26263336932151, 2, 597 -1.06141115981725, 1.28907991440387, 3), 598 ncol = 2) 599 600 y.500 <- log( - coeffqpkwad500[1, ] / p^coeffqpkwad500[2, ] ) 601 y.875 <- log( - coeffqpkwad875[1, ] / p^coeffqpkwad875[2, ] ) 602 603 A.500 <- cbind(1, - log(coeffqpkwad500[3, ] * p^2)) 604 A.875 <- cbind(1, - log(coeffqpkwad875[3, ] * p^2)) 605 coeffic.500 <- solve(A.500, y.500) 606 coeffic.875 <- solve(A.875, y.875) 607 fp.500.n <- 1 - exp(coeffic.500[1]) / n^coeffic.500[2] 608 fp.875.n <- 1 - exp(coeffic.875[1]) / n^coeffic.875[2] 609 } 610 else if(p == 2) { 611 fp.500.n <- 1 - exp( 0.673292623522027) / n^0.691365864961895 612 fp.875.n <- 1 - exp( 0.446537815635445) / n^1.06690782995919 613 } else if(p == 1) { 614 fp.500.n <- 1 - exp( 0.262024211897096) / n^0.604756680630497 615 fp.875.n <- 1 - exp(-0.351584646688712) / n^1.01646567502486 616 } 617 618 ## VT:18.04.2007 - use simulated correction factors for several p and n: 619 ## p in [1, 20] n in [2*p, ...] 620 if(alpha == 0.5 && !is.na(fp.x <- MCDcnp2s$sim.0(p, n))) 621 fp.500.n <- 1/fp.x 622 623 fp.alpha.n <- 624 if(alpha <= 0.875) 625 fp.500.n + (fp.875.n - fp.500.n)/0.375 * (alpha - 0.5) 626 else ## 0.875 < alpha <= 1 627 fp.875.n + (1 - fp.875.n)/0.125 * (alpha - 0.875) 628 629 1/fp.alpha.n 630} ## end{.MCDcnp2 } 631 632MCDcnp2.rew <- # <- *not* exported, but currently used in pkg rrcovNA 633##' Finite sample correction factor for *REW*eighted MCD 634.MCDcnp2.rew <- function(p, n, alpha) 635{ 636 stopifnot(0 <= alpha, alpha <= 1, length(alpha) == 1) 637 638 if(p > 2) { 639 ## "alfaq" "betaq" "qwaarden" 640 coeffrewqpkwad875 <- matrix(c(-0.544482443573914, 1.25994483222292, 2, 641 -0.343791072183285, 1.25159004257133, 3), 642 ncol = 2) 643 coeffrewqpkwad500 <- matrix(c(-1.02842572724793, 1.67659883081926, 2, 644 -0.26800273450853, 1.35968562893582, 3), 645 ncol = 2) 646 647 y.500 <- log( - coeffrewqpkwad500[1, ] / p^ coeffrewqpkwad500[2, ] ) 648 y.875 <- log( - coeffrewqpkwad875[1, ] / p^ coeffrewqpkwad875[2, ] ) 649 650 A.500 <- cbind(1, - log(coeffrewqpkwad500[3, ] * p^2)) 651 coeffic.500 <- solve(A.500, y.500) 652 A.875 <- cbind(1, - log(coeffrewqpkwad875[3, ] * p^2)) 653 coeffic.875 <- solve(A.875, y.875) 654 fp.500.n <- 1 - exp(coeffic.500[1]) / n^ coeffic.500[2] 655 fp.875.n <- 1 - exp(coeffic.875[1]) / n^ coeffic.875[2] 656 } 657 else if(p == 2) { 658 fp.500.n <- 1 - exp( 3.11101712909049 ) / n^ 1.91401056721863 659 fp.875.n <- 1 - exp( 0.79473550581058 ) / n^ 1.10081930350091 660 } else if(p == 1) { 661 fp.500.n <- 1 - exp( 1.11098143415027 ) / n^ 1.5182890270453 662 fp.875.n <- 1 - exp( -0.66046776772861) / n^ 0.88939595831888 663 } 664 665 ## VT:18.04.2007 - use simulated correction factors for several p and n: 666 ## p in [1, 20] n in [2*p, ...] 667 if(alpha == 0.5 && !is.na(fp.x <- MCDcnp2s$sim.rew(p, n))) 668 fp.500.n <- 1/fp.x 669 670 fp.alpha.n <- 671 if(alpha <= 0.875) 672 fp.500.n + (fp.875.n - fp.500.n)/0.375 * (alpha - 0.5) 673 else ## 0.875 < alpha <= 1 674 fp.875.n + (1 - fp.875.n)/0.125 * (alpha - 0.875) 675 676 1/fp.alpha.n 677} ## end{.MCDcnp2.rew } 678 679 680.fastmcd <- function(x, h, nsamp, nmini, kmini, trace = 0) 681{ 682 dx <- dim(x) 683 n <- dx[1] 684 p <- dx[2] 685 686 ## parameters for partitioning {equal to those in Fortran !!} 687 ## kmini <- 5 688 ## nmini <- 300 689 stopifnot(length(kmini <- as.integer(kmini)) == 1, kmini >= 2L, 690 is.1num(nmini), is.finite(nmaxi <- as.double(nmini)*kmini), 691 nmaxi * p < .Machine$integer.max) 692 nmaxi <- as.integer(nmaxi) 693 km10 <- 10*kmini 694 695 ## vt::03.02.2006 - added options "best" and "exact" for nsamp 696 ## 697 nLarge <- 100000 # was 5000 before Nov.2009 -- keep forever now; user can say "exact" 698 if(is.numeric(nsamp) && (nsamp < 0 || nsamp == 0 && p > 1)) { 699 nsamp <- -1 700 } else if(nsamp == "exact" || nsamp == "best") { 701 if(n > 2*nmini-1) { 702 warning("Options 'best' and 'exact' not allowed for n greater than 2*nmini-1 =", 703 2*nmini-1,".\nUsing default.\n") 704 nsamp <- -1 705 } else { 706 myk <- p + 1 ## was 'p'; but p+1 ("nsel = nvar+1") is correct 707 nall <- choose(n, myk) 708 msg <- paste("subsets of size", myk, "out of", n) 709 if(nall > nLarge && nsamp == "best") { 710 nsamp <- nLarge 711 warning("'nsamp = \"best\"' allows maximally ", 712 format(nLarge, scientific=FALSE), 713 " subsets;\ncomputing these ", msg, 714 immediate. = TRUE) 715 } else { ## "exact" or ("best" & nall < nLarge) 716 nsamp <- 0 ## all subsamples -> special treatment in Fortran 717 if(nall > nLarge) { 718 msg <- paste("Computing all", nall, msg) 719 if(nall > 10*nLarge) 720 warning(msg, "\n This may take a", 721 if(nall/nLarge > 100) " very", " long time!\n", 722 immediate. = TRUE) 723 else message(msg) 724 } 725 } 726 } 727 } 728 729 if(!is.numeric(nsamp) || nsamp == -1) { # still not defined 730 ## set it to the default : 731 nsamp.def <- rrcov.control()$nsamp 732 warning(gettextf("Invalid number of trials nsamp=%s. Using default nsamp=%d.", 733 format(nsamp), nsamp.def), 734 domain=NA) 735 nsamp <- nsamp.def 736 } 737 738 if(nsamp > (mx <- .Machine$integer.max)) { 739 warning("nsamp > i_max := maximal integer -- not allowed;\n", 740 " set to i_max = ", mx) 741 nsamp <- mx 742 } 743 744 ## Allocate temporary storage for the Fortran implementation, 745 ## directly in the .Fortran() call. 746 ## (if we used C + .Call() we would allocate all there, and be quite faster!) 747 748 .Fortran(rffastmcd, 749 x = if(is.double(x)) x else as.double(x), 750 n = as.integer(n), 751 p = as.integer(p), ## = 'nvar' in Fortran 752 nhalff = as.integer(h), 753 nsamp = as.integer(nsamp), # = 'krep' 754 nmini = as.integer(nmini), 755 kmini = kmini, 756 initcovariance = double(p * p), 757 initmean = double(p), 758 best = rep.int(as.integer(10000), h), 759 mcdestimate = double(1), ## = 'det' 760 weights = integer(n), 761 exactfit = integer(1), # output indicator: 0: ok; 1: ..., 2: .. 762 coeff = matrix(double(5 * p), nrow = 5, ncol = p), ## plane 763 kount = integer(1), 764 adjustcov = double(p * p), ## used in ltsReg() ! 765 ## integer(1), ## << 'seed' no longer used 766 temp = integer(n), 767 index1 = integer(n), 768 index2 = integer(n), 769 indexx = integer(n), 770 nmahad = double(n), 771 ndist = double(n), 772 am = double(n), 773 am2 = double(n), 774 slutn = double(n), 775 776 med = double(p), 777 mad = double(p), 778 sd = double(p), 779 means = double(p), 780 bmeans= double(p), 781 w = double(p), 782 fv1 = double(p), 783 fv2 = double(p), 784 785 rec = double(p+1), 786 sscp1 = double((p+1)*(p+1)), 787 cova1 = double(p * p), 788 corr1 = double(p * p), 789 cinv1 = double(p * p), 790 cova2 = double(p * p), 791 cinv2 = double(p * p), 792 z = double(p * p), 793 794 cstock = double(10 * p * p), # (10,nvmax2) 795 mstock = double(10 * p), # (10,nvmax) 796 c1stock = double(km10 * p * p), # (km10,nvmax2) 797 m1stock = double(km10 * p), # (km10,nvmax) 798 dath = double(nmaxi * p), # (nmaxi,nvmax) 799 800 cutoff = qchisq(0.975, p), 801 chimed = qchisq(0.5, p), 802 i.trace= as.integer(trace) 803 )[ ## keep the following ones: 804 c("initcovariance", "initmean", "best", "mcdestimate", 805 "weights", "exactfit", "coeff", "kount", "adjustcov") ] 806} 807 808## 809## VT:18.04.2007 - use simulated correction factors for several p and n 810## and alpha = 1/2 (the default in rrcov.control()) 811## ~~~~~~~~~~~ 812## p in [1, 20] n in [2*p, ...] 813## see the modifications in.MCDcnp2() and.MCDcnp2.rew 814## 815 816## VT::08.06.2007 - fixed the simulated values (especially for p=1) 817## VT::11.05.2007 - reduce the usage of the simulated correction factors to only those that 818## are definitvely wrong (negative or very large). This is done by: 819## a) reducing p.max 820## b) reducing n.max 821## NB: In general, "wrong" are the factors for the reweighted matrix, but whenever a simulated 822## value for the reweighted is used, the corresponding simulated must be used for the raw too. 823## 824 825## MM::2014-04 : 826MCDcnp2s <- local({ 827 p.min <- 1L 828 p.max <- 9L # was 20 829 ncol <- 20L # the number of column in the matrices 830 n.min <- as.integer( 831### p = 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 832 c(1, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40)) 833 n.max <- as.integer( 834 c(2, 6, 10, 13, 16, 18, 20, 20, 20, 40, 42, 44, 46, 48, 50, 52, 54, 56, 58, 60)) 835##was c(22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 54, 56, 58, 60) 836## these are the right (simulated) values for n.max 837 838 n.min.rew <- n.min 839 n.max.rew <- n.max 840 841 m.0 <- matrix( 842 c(1, 3.075819, 1.515999, 2.156169, 1.480742, 1.765485, 1.460206, 1.603707, 1.427429, 1.504712, 1.334528, 1.48297, 1.355308, 1.383867, 1.319241, 1.36065, 1.307467, 1.365596, 1.255259, 1.352741, 1.239381, 3.15342, 1.799889, 2.258497, 1.688312, 1.906779, 1.548203, 1.724785, 1.500873, 1.573442, 1.417137, 1.540805, 1.395945, 1.472596, 1.394247, 1.377487, 1.337394, 1.369354, 1.333378, 1.3181, 1.313813, 1.315528, 2.12777, 2.718898, 1.993509, 2.220433, 1.820585, 1.97782, 1.672455, 1.770151, 1.587478, 1.685352, 1.539295, 1.584536, 1.499487, 1.50702, 1.41952, 1.449058, 1.393042, 1.432999, 1.369964, 1.400997, 1.333824, 2.950549, 2.145387, 2.382224, 1.927077, 2.032489, 1.8371, 1.877833, 1.710891, 1.756053, 1.620778, 1.657761, 1.558978, 1.56257, 1.508633, 1.534406, 1.46709, 1.468734, 1.432529, 1.455283, 1.386975, 1.417532, 2.229573, 2.494447, 2.016117, 2.190061, 1.877996, 1.978964, 1.767284, 1.836948, 1.677372, 1.743316, 1.616383, 1.655964, 1.55484, 1.594831, 1.502185, 1.543723, 1.467005, 1.491123, 1.44402, 1.446915, 1.401578, 2.580264, 2.109121, 2.240741, 1.944719, 2.043397, 1.821808, 1.89725, 1.748788, 1.786988, 1.659333, 1.697012, 1.610622, 1.616503, 1.538529, 1.562024, 1.499964, 1.529344, 1.474519, 1.483264, 1.441552, 1.434448, 2.165233, 2.320281, 2.007836, 2.086471, 1.884052, 1.950563, 1.76926, 1.843328, 1.708941, 1.741039, 1.627206, 1.644755, 1.580563, 1.593402, 1.527312, 1.568418, 1.501462, 1.502542, 1.464583, 1.467921, 1.431141, 2.340443, 2.048262, 2.161097, 1.926082, 1.995422, 1.81446, 1.853165, 1.738533, 1.784456, 1.679444, 1.696463, 1.612931, 1.629483, 1.548186, 1.580026, 1.52198, 1.531111, 1.482914, 1.484824, 1.442726, 1.447838, 2.093386, 2.185793, 1.948989, 2.02804, 1.867137, 1.907732, 1.771923, 1.800413, 1.691612, 1.720603, 1.642705, 1.649769, 1.589028, 1.598955, 1.539759, 1.55096, 1.503965, 1.50703, 1.471349, 1.469791, 1.436959, 2.218315, 1.997369, 2.041128, 1.887059, 1.928524, 1.79626, 1.827538, 1.716748, 1.735696, 1.658329, 1.664211, 1.599286, 1.611511, 1.553925, 1.562637, 1.516805, 1.529894, 1.476064, 1.482474, 1.453253, 1.458467, 2.0247, 2.07899, 1.921976, 1.949376, 1.824629, 1.851671, 1.744713, 1.765647, 1.683525, 1.685592, 1.625113, 1.624961, 1.571921, 1.581223, 1.535257, 1.537464, 1.497165, 1.504879, 1.468682, 1.469319, 1.448344, 2.092315, 1.941412, 1.969843, 1.844093, 1.866133, 1.766145, 1.783829, 1.703613, 1.709714, 1.646078, 1.654264, 1.594523, 1.598488, 1.545105, 1.555356, 1.514627, 1.521353, 1.483958, 1.487677, 1.449191, 1.459721, 1.958987, 1.985144, 1.87739, 1.879643, 1.786823, 1.799642, 1.720015, 1.724688, 1.663539, 1.662997, 1.609267, 1.615124, 1.56746, 1.562026, 1.520586, 1.52503, 1.493008, 1.502496, 1.471983, 1.468546, 1.435064, 1.994706, 1.880348, 1.894254, 1.805827, 1.815965, 1.744296, 1.743389, 1.665481, 1.681644, 1.624466, 1.626109, 1.584028, 1.5818, 1.54376, 1.547237, 1.504878, 1.515087, 1.479032, 1.47936, 1.450758, 1.45073, 1.892685, 1.91087, 1.825301, 1.827176, 1.745363, 1.746115, 1.693373, 1.701692, 1.648247, 1.637112, 1.594648, 1.592013, 1.554849, 1.55013, 1.522186, 1.520901, 1.492606, 1.493072, 1.460868, 1.46733, 1.440956, 1.92771, 1.835696, 1.841979, 1.775991, 1.766092, 1.703807, 1.708791, 1.654985, 1.655917, 1.602388, 1.611867, 1.570765, 1.573368, 1.53419, 1.529033, 1.506767, 1.503596, 1.481126, 1.471806, 1.444917, 1.451682, 1.850262, 1.855034, 1.778997, 1.789995, 1.718871, 1.717326, 1.667357, 1.666291, 1.619743, 1.631475, 1.582624, 1.58766, 1.546302, 1.545063, 1.512222, 1.517888, 1.489127, 1.487271, 1.466722, 1.463618, 1.444137, 1.8709, 1.794033, 1.80121, 1.736376, 1.740201, 1.673776, 1.682541, 1.638153, 1.642294, 1.604417, 1.597721, 1.559534, 1.559108, 1.533942, 1.529348, 1.499517, 1.501586, 1.473147, 1.473031, 1.457615, 1.452348, 1.805753, 1.812952, 1.746549, 1.747222, 1.696924, 1.694957, 1.652157, 1.650568, 1.607807, 1.613666, 1.577295, 1.570712, 1.543704, 1.538272, 1.515369, 1.517113, 1.487451, 1.491593, 1.464514, 1.464658, 1.439359, 1.823222, 1.758781, 1.767358, 1.70872, 1.712926, 1.666956, 1.667838, 1.62077, 1.621445, 1.592891, 1.58549, 1.55603, 1.559042, 1.521501, 1.523342, 1.499913, 1.501937, 1.473359, 1.472522, 1.452613, 1.452448), 843 ncol = ncol) 844 845 m.rew <- matrix( 846 c(1, 0.984724, 0.970109, 0.978037, 0.979202, 0.982933, 1.001461, 1.026651, 0.981233, 1.011895, 1.017499, 0.964323, 1.026574, 1.006594, 0.980194, 1.009828, 0.998083, 0.966173, 1.009942, 0.99916, 1.021521, 2.216302, 1.418526, 1.635601, 1.31402, 1.33975, 1.251798, 1.210917, 1.133114, 1.150666, 1.138732, 1.096822, 1.076489, 1.058343, 1.045746, 1.036743, 1.008929, 1.049537, 1.028148, 1.027297, 1.020578, 1.00074, 1.73511, 2.06681, 1.545905, 1.659655, 1.456835, 1.47809, 1.331966, 1.334229, 1.231218, 1.220443, 1.198143, 1.193965, 1.142156, 1.146231, 1.124661, 1.112719, 1.089973, 1.070606, 1.082681, 1.061243, 1.053191, 2.388892, 1.847626, 1.96998, 1.630723, 1.701272, 1.521008, 1.553057, 1.382168, 1.414555, 1.326982, 1.321403, 1.265207, 1.264856, 1.200418, 1.21152, 1.17531, 1.168536, 1.140586, 1.14457, 1.111392, 1.112031, 1.968153, 2.168931, 1.784373, 1.894409, 1.667912, 1.693007, 1.545176, 1.582428, 1.45319, 1.480559, 1.371611, 1.358541, 1.330235, 1.30264, 1.257518, 1.244156, 1.221907, 1.22455, 1.178965, 1.177855, 1.166319, 2.275891, 1.866587, 2.014249, 1.750567, 1.829363, 1.650019, 1.689043, 1.562539, 1.561359, 1.473378, 1.488554, 1.411097, 1.416527, 1.35117, 1.361044, 1.30205, 1.299037, 1.250265, 1.260083, 1.218665, 1.236027, 1.95771, 2.074066, 1.847385, 1.905408, 1.71393, 1.768425, 1.63908, 1.67234, 1.564992, 1.562337, 1.49229, 1.499573, 1.420813, 1.424067, 1.383947, 1.378726, 1.33062, 1.330071, 1.279404, 1.295302, 1.263947, 2.164121, 1.871024, 1.979485, 1.782417, 1.84489, 1.706023, 1.734857, 1.622782, 1.634869, 1.55196, 1.554423, 1.482325, 1.509195, 1.440726, 1.436328, 1.386335, 1.396277, 1.347939, 1.346732, 1.310242, 1.309371, 1.938822, 2.050409, 1.834863, 1.882536, 1.737494, 1.761608, 1.65742, 1.687579, 1.591863, 1.60158, 1.520982, 1.535234, 1.470649, 1.486485, 1.42892, 1.435574, 1.384132, 1.382329, 1.343281, 1.346581, 1.315111, 2.063894, 1.880094, 1.907246, 1.78278, 1.806648, 1.6952, 1.720922, 1.63084, 1.635274, 1.565423, 1.56171, 1.512015, 1.4986, 1.463903, 1.456588, 1.422856, 1.407325, 1.376724, 1.373923, 1.346464, 1.34259, 1.898389, 1.950406, 1.812053, 1.849175, 1.72649, 1.737651, 1.646719, 1.655112, 1.587601, 1.597894, 1.539877, 1.53329, 1.495054, 1.490548, 1.445249, 1.446037, 1.410272, 1.412274, 1.375797, 1.369604, 1.341232, 1.992488, 1.830452, 1.857314, 1.758686, 1.763822, 1.683215, 1.679543, 1.619269, 1.608512, 1.565, 1.562282, 1.498869, 1.51325, 1.470912, 1.464654, 1.427573, 1.439301, 1.402308, 1.391006, 1.37074, 1.367573, 1.855502, 1.891242, 1.77513, 1.790618, 1.706443, 1.713098, 1.642896, 1.636577, 1.580366, 1.581752, 1.542937, 1.531668, 1.487894, 1.492039, 1.460304, 1.449762, 1.4219, 1.420953, 1.390137, 1.388677, 1.360506, 1.908277, 1.802091, 1.806128, 1.723757, 1.727249, 1.659883, 1.670056, 1.605209, 1.611481, 1.558846, 1.551762, 1.512951, 1.511515, 1.468948, 1.476073, 1.441508, 1.434997, 1.412687, 1.406782, 1.380452, 1.375924, 1.811415, 1.822311, 1.740544, 1.739355, 1.68127, 1.685342, 1.620281, 1.622572, 1.579611, 1.570103, 1.529881, 1.530097, 1.490041, 1.4947, 1.457329, 1.456344, 1.423363, 1.428653, 1.399988, 1.390069, 1.376594, 1.837723, 1.76039, 1.771031, 1.697404, 1.690915, 1.634409, 1.63713, 1.589594, 1.586521, 1.552974, 1.545571, 1.505923, 1.512794, 1.477833, 1.477821, 1.444241, 1.44452, 1.419258, 1.421297, 1.394924, 1.389393, 1.779716, 1.781271, 1.706031, 1.71224, 1.655099, 1.654284, 1.608878, 1.605955, 1.565683, 1.565938, 1.523594, 1.531235, 1.492749, 1.486786, 1.457635, 1.461416, 1.432472, 1.430164, 1.404441, 1.400021, 1.378273, 1.798932, 1.735577, 1.727031, 1.671049, 1.677601, 1.624427, 1.617626, 1.579533, 1.579987, 1.544635, 1.538715, 1.504538, 1.50726, 1.477163, 1.477084, 1.450861, 1.444496, 1.428416, 1.422813, 1.400185, 1.39552, 1.750193, 1.752145, 1.690365, 1.692051, 1.642391, 1.63858, 1.600144, 1.596401, 1.558305, 1.555932, 1.525968, 1.522984, 1.491563, 1.492554, 1.467575, 1.45786, 1.437545, 1.430893, 1.413983, 1.409386, 1.391943, 1.762922, 1.701346, 1.704996, 1.6556, 1.655548, 1.611964, 1.615219, 1.569103, 1.571079, 1.540617, 1.541602, 1.503791, 1.50195, 1.478069, 1.47678, 1.452458, 1.451732, 1.429144, 1.426547, 1.40363, 1.402647), 847 ncol = ncol) 848 849 rm(ncol) 850 list( 851 sim.0 = function(p, n) { 852 p. <- p - p.min + 1L 853 if(p.min <= p && p <= p.max && 854 n.min[p.] <= n && n <= n.max[p.]) { 855 nind <- n - n.min[p.] + 1L 856 m.0[nind, p.] 857 ##= 858 } else NA 859 }, 860 sim.rew = function(p, n) { 861 p. <- p - p.min + 1L 862 if(p.min <= p && p <= p.max && 863 n.min.rew[p.] <= n && n <= n.max.rew[p.]) { 864 nind <- n - n.min.rew[p.] + 1L 865 m.rew[nind, p.] 866 ##=== 867 } else NA 868 }) 869}) ## end{MCDcnp2s} 870 871if(FALSE) { ## For experimentation: 872 ls.str( ee <- environment(MCDcnp2s$sim.0) ) 873 matplot(1:21, ee$m.0, type = "o", xlab = "p - p.min + 1") 874} 875