@@ -140,8 +140,10 @@ plot.mclustBootstrapLRT <- function(x, G = 1, hist.col = "grey", hist.border = "
140140# Bootstrap inference (standard errors and percentile confidence intervals)
141141#
142142
143- MclustBootstrap <- function (object , nboot = 999 , type = c(" bs" " wlbs" " pb" " jk"
144- max.nonfit = 10 * nboot , verbose = interactive(), ... )
143+ MclustBootstrap <- function (object , nboot = 999 ,
144+ type = c(" bs" " wlbs" " pb" " jk"
145+ alpha = 1 , max.nonfit = 10 * nboot ,
146+ verbose = interactive(), ... )
145147{
146148
147149 stopifnot(" object must be of class 'Mclust' or 'densityMclust'" =
@@ -195,7 +197,8 @@ MclustBootstrap <- function(object, nboot = 999, type = c("bs", "wlbs", "pb", "j
195197 },
196198 " wlbs" =
197199 {
198- w <- rexp(n )
200+ # w ~ Dirichlet(alpha, ..., alpha)
201+ w <- rgamma(n , shape = as.numeric(alpha ), rate = 1 )
199202 # w <- w/mean(w)
200203 w <- w / max(w )
201204 mod.boot <- try(do.call(" me.weighted"
@@ -279,60 +282,64 @@ summary.MclustBootstrap <- function(object, what = c("se", "ci", "ave"), conf.le
279282 G <- dims [3 ]
280283
281284 switch (what ,
282- " se" = { out <- list (pro = apply(object $ pro , 2 , sd , na.rm = TRUE ),
283- mean = apply(object $ mean , c(2 ,3 ), sd , na.rm = TRUE ),
284- variance = apply(object $ variance , c(2 ,3 ,4 ), sd , na.rm = TRUE ))
285- if (object $ type == " jk"
286- out <- lapply(out , function (x )
287- sqrt(x ^ 2 * (nboot - object $ nonfit - 1 )^ 2 / (nboot - object $ nonfit )))
288- },
289- " ave" = { out <- list (pro = apply(object $ pro , 2 , mean , na.rm = TRUE ),
290- mean = apply(object $ mean , c(2 ,3 ), mean , na.rm = TRUE ),
291- variance = apply(object $ variance , c(2 ,3 ,4 ), mean , na.rm = TRUE ))
292- },
293- " ci" = { levels <- c((1 - conf.level )/ 2 , (1 + conf.level )/ 2 )
294- if (object $ type == " jk"
295- { # bias-corrected ci based on normal-approximation
296- ave <- list (pro = apply(object $ pro , 2 , mean , na.rm = TRUE ),
297- mean = apply(object $ mean , c(2 ,3 ), mean , na.rm = TRUE ),
298- variance = t(sapply(seq.int(d ), function (j )
299- apply(object $ variance [,j ,j ,], 2 , mean , na.rm = TRUE ),
300- simplify = " array"
301- se <- list (pro = apply(object $ pro , 2 , sd , na.rm = TRUE ),
302- mean = apply(object $ mean , c(2 ,3 ), sd , na.rm = TRUE ),
303- variance = t(sapply(seq.int(d ), function (j )
304- apply(object $ variance [,j ,j ,], 2 , sd , na.rm = TRUE ),
305- simplify = " array"
306- se <- lapply(se , function (x )
307- sqrt(x ^ 2 * (nboot - object $ nonfit - 1 )^ 2 / (nboot - object $ nonfit )))
308- zq <- qnorm(max(levels ))
309- lnames <- paste0(formatC(levels * 100 , format = " fg" width = 1 ,
310- digits = getOption(" digits" " %"
311- # the code above mimic stats:::format_perc(levels) which can't be used
312- # because format_perc is not exported from stats
313- out <- list (pro = array (as.double(NA ), c(2 ,G ),
314- dimnames = list (lnames , 1 : G )),
315- mean = array (as.double(NA ), dim = c(2 ,d ,G ),
316- dimnames = list (lnames , 1 : d , 1 : G )),
317- variance = array (as.double(NA ), dim = c(2 ,d ,G ),
318- dimnames = list (lnames , 1 : d , 1 : G )))
319- out $ pro [1 ,] <- ave $ pro - zq * se $ pro
320- out $ pro [2 ,] <- ave $ pro + zq * se $ pro
321- out $ mean [1 ,,] <- ave $ mean - zq * se $ mean
322- out $ mean [2 ,,] <- ave $ mean + zq * se $ mean
323- out $ variance [1 ,,] <- ave $ variance - zq * se $ variance
324- out $ variance [2 ,,] <- ave $ variance + zq * se $ variance
325- } else
326- { # percentile-based ci
327- out <- list (pro = apply(object $ pro , 2 , quantile , probs = levels , na.rm = TRUE ),
328- mean = apply(object $ mean , c(2 ,3 ), quantile , probs = levels , na.rm = TRUE ))
329- v <- array (as.double(NA ), dim = c(2 ,d ,G ),
330- dimnames = dimnames(out $ mean ))
331- for (j in seq.int(d ))
332- v [,j ,] <- apply(object $ variance [,j ,j ,], 2 , quantile , probs = levels , na.rm = TRUE )
333- out $ variance <- v
334- }
335- }
285+ " se" = {
286+ out <- list (pro = apply(object $ pro , 2 , sd , na.rm = TRUE ),
287+ mean = apply(object $ mean , c(2 ,3 ), sd , na.rm = TRUE ),
288+ variance = apply(object $ variance , c(2 ,3 ,4 ), sd , na.rm = TRUE ))
289+ if (object $ type == " jk"
290+ out <- lapply(out , function (x )
291+ sqrt(x ^ 2 * (nboot - object $ nonfit - 1 )^ 2 / (nboot - object $ nonfit )))
292+ },
293+ " ave" = {
294+ out <- list (pro = apply(object $ pro , 2 , mean , na.rm = TRUE ),
295+ mean = apply(object $ mean , c(2 ,3 ), mean , na.rm = TRUE ),
296+ variance = apply(object $ variance , c(2 ,3 ,4 ), mean , na.rm = TRUE ))
297+ },
298+ " ci" = {
299+ levels <- c((1 - conf.level )/ 2 , (1 + conf.level )/ 2 )
300+ if (object $ type == " jk"
301+ { # bias-corrected ci based on normal-approximation
302+ ave <- list (pro = apply(object $ pro , 2 , mean , na.rm = TRUE ),
303+ mean = apply(object $ mean , c(2 ,3 ), mean , na.rm = TRUE ),
304+ variance = t(sapply(seq.int(d ), function (j )
305+ apply(object $ variance [,j ,j ,], 2 , mean , na.rm = TRUE ),
306+ simplify = " array"
307+ se <- list (pro = apply(object $ pro , 2 , sd , na.rm = TRUE ),
308+ mean = apply(object $ mean , c(2 ,3 ), sd , na.rm = TRUE ),
309+ variance = t(sapply(seq.int(d ), function (j )
310+ apply(object $ variance [,j ,j ,], 2 , sd , na.rm = TRUE ),
311+ simplify = " array"
312+ se <- lapply(se , function (x )
313+ sqrt(x ^ 2 * (nboot - object $ nonfit - 1 )^ 2 / (nboot - object $ nonfit )))
314+ zq <- qnorm(max(levels ))
315+ lnames <- paste0(formatC(levels * 100 , format = " fg" width = 1 ,
316+ digits = getOption(" digits" " %"
317+ # the code above mimic stats:::format_perc(levels) which can't be used
318+ # because format_perc is not exported from stats
319+ out <- list (pro = array (as.double(NA ), c(2 ,G ),
320+ dimnames = list (lnames , 1 : G )),
321+ mean = array (as.double(NA ), dim = c(2 ,d ,G ),
322+ dimnames = list (lnames , 1 : d , 1 : G )),
323+ variance = array (as.double(NA ), dim = c(2 ,d ,G ),
324+ dimnames = list (lnames , 1 : d , 1 : G )))
325+ out $ pro [1 ,] <- ave $ pro - zq * se $ pro
326+ out $ pro [2 ,] <- ave $ pro + zq * se $ pro
327+ out $ mean [1 ,,] <- ave $ mean - zq * se $ mean
328+ out $ mean [2 ,,] <- ave $ mean + zq * se $ mean
329+ out $ variance [1 ,,] <- ave $ variance - zq * se $ variance
330+ out $ variance [2 ,,] <- ave $ variance + zq * se $ variance
331+ } else
332+ { # percentile-based ci
333+ out <- list (pro = apply(object $ pro , 2 , quantile , probs = levels , na.rm = TRUE ),
334+ mean = apply(object $ mean , c(2 ,3 ), quantile , probs = levels , na.rm = TRUE ))
335+ v <- array (as.double(NA ), dim = c(2 ,d ,G ),
336+ dimnames = dimnames(out $ mean ))
337+ for (j in seq.int(d ))
338+ v [,j ,] <- apply(object $ variance [,j ,j ,,drop = FALSE ], 4 ,
339+ quantile , probs = levels , na.rm = TRUE )
340+ out $ variance <- v
341+ }
342+ }
336343 )
337344
338345 obj <- append(object [c(" modelName" " G" " nboot" " type"
@@ -399,17 +406,18 @@ print.summary.MclustBootstrap <- function(x, digits = getOption("digits"), ...)
399406 invisible (x )
400407}
401408
402- plot.MclustBootstrap <- function (x , what = c(" pro" " mean" " var" show.parest = TRUE , show.confint = TRUE , hist.col = " grey" hist.border = " lightgrey" breaks = " Sturges " col = " forestgreen" lwd = 2 , lty = 3 , xlab = NULL , xlim = NULL , ylim = NULL , ... )
409+ plot.MclustBootstrap <- function (x , what = c(" pro" " mean" " var" show.parest = TRUE , show.confint = TRUE , hist.col = " grey" hist.border = " lightgrey" breaks = NA , col = " forestgreen" lwd = 2 , lty = 3 , xlab = NULL , xlim = NULL , ylim = NULL , ... )
403410{
404411 object <- x # Argh. Really want to use object anyway
405412 what <- match.arg(what , choices = eval(formals(plot.MclustBootstrap )$ what ))
406413 par <- object $ parameters
407414 d <- dim(object $ mean )[2 ]
408415 varnames <- rownames(par $ mean )
409416 if (show.confint )
410- { ci <- summary(object , what = " ci" ... )
411- ave <- summary(object , what = " ave" ... )
412- }
417+ {
418+ ci <- summary(object , what = " ci" ... )
419+ ave <- summary(object , what = " ave" ... )
420+ }
413421
414422 histBoot <- function (boot , stat , ci , ave , breaks , xlim , ylim , xlab , ... )
415423 {
@@ -426,48 +434,54 @@ plot.MclustBootstrap <- function(x, what = c("pro", "mean", "var"), show.parest
426434 }
427435
428436 switch (what ,
429- " pro" = { xlim <- range(if (is.null(xlim )) pretty(object $ pro ) else xlim )
430- for (k in 1 : object $ G )
431- histBoot(object $ pro [,k ], breaks = breaks ,
432- stat = par $ pro [k ],
433- ci = ci $ pro [,k ],
434- ave = ave $ pro [k ],
435- xlim = xlim , ylim = ylim ,
436- xlab = ifelse(is.null(xlab ),
437- paste(" Mix. prop. for comp." k ),
438- xlab ))
437+ " pro" = {
438+ xlim <- range(if (is.null(xlim )) pretty(object $ pro ) else xlim )
439+ for (k in 1 : object $ G )
440+ histBoot(object $ pro [,k ],
441+ breaks = if (is.na(breaks )) nclass.numpy else breaks ,
442+ stat = par $ pro [k ],
443+ ci = ci $ pro [,k ],
444+ ave = ave $ pro [k ],
445+ xlim = xlim , ylim = ylim ,
446+ xlab = ifelse(is.null(xlab ),
447+ paste(" Mix. prop. for comp." k ),
448+ xlab ))
439449 },
440- " mean" = { isNull_xlim <- is.null(xlim )
441- for (j in 1 : d )
442- { xlim <- range(if (isNull_xlim ) pretty(object $ mean [,j ,])
443- else xlim )
444- for (k in 1 : object $ G )
445- histBoot(object $ mean [,j ,k ], breaks = breaks ,
446- stat = par $ mean [j ,k ],
447- ci = ci $ mean [,j ,k ],
448- ave = ave $ mean [j ,k ],
449- xlim = xlim , ylim = ylim ,
450- xlab = ifelse(is.null(xlab ),
451- paste(varnames [j ], " mean for comp." k ),
452- xlab ))
453- }
450+ " mean" = {
451+ isNull_xlim <- is.null(xlim )
452+ for (j in 1 : d )
453+ {
454+ xlim <- range(if (isNull_xlim ) pretty(object $ mean [,j ,]) else xlim )
455+ for (k in 1 : object $ G )
456+ histBoot(object $ mean [,j ,k ],
457+ breaks = if (is.na(breaks )) nclass.numpy else breaks ,
458+ stat = par $ mean [j ,k ],
459+ ci = ci $ mean [,j ,k ],
460+ ave = ave $ mean [j ,k ],
461+ xlim = xlim , ylim = ylim ,
462+ xlab = ifelse(is.null(xlab ),
463+ paste(varnames [j ], " mean for comp." k ),
464+ xlab ))
465+ }
454466 },
455- " var" = { isNull_xlim <- is.null(xlim )
456- for (j in 1 : d )
457- { xlim <- range(if (isNull_xlim ) pretty(object $ variance [,j ,j ,])
458- else xlim )
459- for (k in 1 : object $ G )
460- histBoot(object $ variance [,j ,j ,k ], breaks = breaks ,
461- stat = par $ variance [j ,j ,k ],
462- ci = ci $ variance [,j ,k ],
463- ave = ave $ variance [j ,k ],
464- xlim = xlim , ylim = ylim ,
465- xlab = ifelse(is.null(xlab ),
466- paste(varnames [j ], " var. for comp." k ),
467- xlab ))
468- }
467+ " var" = {
468+ isNull_xlim <- is.null(xlim )
469+ for (j in 1 : d )
470+ {
471+ xlim <- range(if (isNull_xlim ) pretty(object $ variance [,j ,j ,]) else xlim )
472+ for (k in 1 : object $ G )
473+ histBoot(object $ variance [,j ,j ,k ],
474+ breaks = if (is.na(breaks )) nclass.numpy else breaks ,
475+ stat = par $ variance [j ,j ,k ],
476+ ci = ci $ variance [,j ,k ],
477+ ave = ave $ variance [j ,j ,k ],
478+ xlim = xlim , ylim = ylim ,
479+ xlab = ifelse(is.null(xlab ),
480+ paste(varnames [j ], " var. for comp." k ),
481+ xlab ))
482+ }
469483 }
470- )
484+ )
471485 invisible ()
472486}
473487
0 commit comments