[Archetypes-commits] r25 - in branches/pkg-robust: . R sandbox
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Feb 22 18:20:18 CET 2010
Author: manuel
Date: 2010-02-22 18:20:18 +0100 (Mon, 22 Feb 2010)
New Revision: 25
Added:
branches/pkg-robust/sandbox/
branches/pkg-robust/sandbox/demo.R
Modified:
branches/pkg-robust/R/archetypes-diagplots.R
branches/pkg-robust/R/archetypes-kit-blocks.R
branches/pkg-robust/R/archetypes-kit.R
branches/pkg-robust/R/archetypes-movie.R
branches/pkg-robust/R/archetypes-plot.R
Log:
(1) diagnostics and (2) step-by-step demo.
Modified: branches/pkg-robust/R/archetypes-diagplots.R
===================================================================
--- branches/pkg-robust/R/archetypes-diagplots.R 2010-02-19 17:26:00 UTC (rev 24)
+++ branches/pkg-robust/R/archetypes-diagplots.R 2010-02-22 17:20:18 UTC (rev 25)
@@ -31,7 +31,11 @@
-rss.diagplot <- function(object, sort = FALSE, ...) {
+rss.diagplot <- function(object, ...) {
+ UseMethod('rss.diagplot')
+}
+
+rss.diagplot.archetypes <- function(object, sort = FALSE, ...) {
y <- rss(object, type = 'single')
x <- seq(length = length(y))
@@ -42,10 +46,58 @@
abline(h = 0, lty = 2, col = gray(0.7), ...)
}
+rss.diagplot.repArchetypes <- function(object, ...) {
+ y <- lapply(object, rss, type = 'single')
+ x <- seq(length = length(y[[1]]))
+ ylim <- range(sapply(y, range))
+
+ plot(x, y[[1]], xlab = 'Index', ylab = 'RSS',
+ ylim = ylim, type = 'n', ...)
+ for ( i in seq(along = y) )
+ lines(x, y[[i]], col = i, ...)
+}
+
+
+
weights.diagplot <- function(object, weights.type, ...) {
y <- weights(object, weights.type)
x <- seq(length = length(y))
- plot(x, y, ylim = c(1, 0), xlab = 'Index', ylab = 'Weights', ...)
+ ylab <- sprintf('%s%s', toupper(substring(weights.type, 1, 1)),
+ substring(weights.type, 2))
+
+ plot(x, y, ylim = c(1, 0), xlab = 'Index', ylab = ylab, ...)
}
+
+
+
+reweights.diagplot <- function(object, highlight = NULL,
+ highlight.col = (seq(length(highlight)) + 1), ...) {
+
+ y <- rev(lapply(object$history, function(x) x[[1]]$reweights))
+ x <- seq(along = y[[1]])
+ col <- rep(1, length = length(x))
+ col[highlight] <- highlight.col
+
+ n <- sqrt(length(y))
+
+ par(mfrow = c(ceiling(n), ceiling(n)), mar = c(0, 0, 0, 0))
+ for ( i in seq(along = y) )
+ plot(x, y[[i]], type = 'p', col = col, xlab = 'Index',
+ ylab = 'Reweights', ...)
+}
+
+
+
+i.reweights.diagplot <- function(object, i, lty = 1,
+ col = (seq(length(i)) + 1), ...) {
+ y <- sapply(object$history, function(x) x[[1]]$reweights[i])
+ y <- apply(y, 1, rev)
+
+ matplot(y, type = 'l', lty = lty, col = col, ylab = 'Reweights',
+ xlab = 'Iterations', ...)
+}
+
+
+
Modified: branches/pkg-robust/R/archetypes-kit-blocks.R
===================================================================
--- branches/pkg-robust/R/archetypes-kit-blocks.R 2010-02-19 17:26:00 UTC (rev 24)
+++ branches/pkg-robust/R/archetypes-kit-blocks.R 2010-02-22 17:20:18 UTC (rev 25)
@@ -300,6 +300,9 @@
#'
#'
bisquare0.reweightsfn <- function(resid) {
+ if ( is.null(resid) )
+ return(NULL)
+
resid <- apply(resid, 2, function(.) sum(abs(.)))
resid0 <- resid < sqrt(.Machine$double.eps)
@@ -309,6 +312,9 @@
}
tricube.reweightsfn <- function(resid) {
+ if ( is.null(resid) )
+ return(NULL)
+
resid <- apply(resid, 2, function(.) sum(abs(.)))
ifelse(resid < 1, (1 - resid^3)^3, 0)
}
Modified: branches/pkg-robust/R/archetypes-kit.R
===================================================================
--- branches/pkg-robust/R/archetypes-kit.R 2010-02-19 17:26:00 UTC (rev 24)
+++ branches/pkg-robust/R/archetypes-kit.R 2010-02-22 17:20:18 UTC (rev 25)
@@ -40,7 +40,8 @@
list(archetypes = as.archetypes(t(family$rescalefn(x, family$undummyfn(x, zs))),
k, alphas = t(alphas), betas = t(betas), rss = rss, kappas = kappas,
zas = t(family$rescalefn(x, family$undummyfn(x, zas))),
- residuals = resid, reweights = reweights))
+ residuals = resid, reweights = reweights,
+ family = list(class = family$class)))
}
printIter <- function(i) {
@@ -67,13 +68,13 @@
betas <- init$betas
alphas <- init$alphas
+ zas <- NULL
zs <- x %*% betas
resid <- zs %*% alphas - x
rss <- family$normfn(resid) / n
- zas <- NULL
- reweights <- NULL
+ reweights <- rep(1, n)
kappas <- c(alphas=kappa(alphas), betas=kappa(betas),
zas=-Inf, zs=kappa(zs))
Modified: branches/pkg-robust/R/archetypes-movie.R
===================================================================
--- branches/pkg-robust/R/archetypes-movie.R 2010-02-19 17:26:00 UTC (rev 24)
+++ branches/pkg-robust/R/archetypes-movie.R 2010-02-22 17:20:18 UTC (rev 25)
@@ -17,8 +17,8 @@
#' @export
movieplot <- function(zs, data, show=c('atypes', 'adata'),
ssleep=0, bsleep=0, postfn=function(iter){}, ...) {
-
- steps <- length(zs$history)
+
+ steps <- length(zs$history)
atypesmovie <- ifelse(show[1] == 'atypes', TRUE, FALSE)
Sys.sleep(ssleep)
@@ -32,7 +32,7 @@
plot(adata(a), ...)
postfn(i)
-
+
Sys.sleep(bsleep)
}
}
@@ -82,7 +82,7 @@
plot(a, data, ...)
Sys.sleep(bsleep)
}
-
+
plot(a, data, ...)
}
@@ -113,7 +113,7 @@
pcplot(a, data, ...)
else
pcplot(adata(a), rx=rx, ...)
-
+
Sys.sleep(bsleep)
}
}
Modified: branches/pkg-robust/R/archetypes-plot.R
===================================================================
--- branches/pkg-robust/R/archetypes-plot.R 2010-02-19 17:26:00 UTC (rev 24)
+++ branches/pkg-robust/R/archetypes-plot.R 2010-02-22 17:20:18 UTC (rev 25)
@@ -86,22 +86,21 @@
plot.weightedArchetypes <- function(x, y,
adata.show = FALSE, data.col = 1, data.pch = 21,
data.bg = gray, link.col.show = TRUE,
+ link.col = data.col, link.lty = 1,
weights.type = 'weights', ...) {
- if ( adata.show ) {
- w <- 1 - weights(x, type = weights.type)
+ if ( !adata.show )
+ return(plot.archetypes(x, y, ...))
- if ( link.col.show ) {
- link.col <- ifelse(w == 1, 1, data.bg(w))
- link.lty <- ifelse(w == 1, 2, 1)
- }
+ w <- 1 - weights(x, type = weights.type)
- plot.archetypes(x, y, adata.show = TRUE, data.pch = data.pch,
- data.col = data.col, data.bg = data.bg(w),
- link.col = link.col, link.lty = link.lty)
+ if ( link.col.show ) {
+ link.col <- ifelse(w == 1, 1, data.bg(w))
+ link.lty <- ifelse(w == 1, 2, 1)
}
- else {
- plot.archetypes(x, y, ...)
- }
+
+ plot.archetypes(x, y, adata.show = TRUE, data.pch = data.pch,
+ data.col = data.col, data.bg = data.bg(w),
+ link.col = link.col, link.lty = link.lty, ...)
}
Added: branches/pkg-robust/sandbox/demo.R
===================================================================
--- branches/pkg-robust/sandbox/demo.R (rev 0)
+++ branches/pkg-robust/sandbox/demo.R 2010-02-22 17:20:18 UTC (rev 25)
@@ -0,0 +1,156 @@
+
+library(archetypes.dev.robust)
+
+
+
+### Archetypes with one outlier; breakdown point simulations ########
+
+data(toy)
+
+
+### Move one data point off the data cloud:
+
+plot(toy, pch = 21, col = 1, bg = c(rep(0, 183), 2, rep(0, 66)))
+
+out.start <- toy[184, ]
+out.end <- c(0, 30)
+
+p184 <- list()
+a184 <- list()
+
+set.seed(1234)
+for ( i in 1:10 ) {
+ p184[[i]] <- out.end - (i/10) * (out.end - out.start)
+ toy[184, ] <- p184[[i]]
+
+ a184[[i]] <- archetypes(toy, 3)
+}
+
+save(p184, a184, file = 'a184.RData')
+load(file = 'a184.RData')
+
+
+plot(rbind(toy, do.call(rbind, p184)), pch = 21, col = 1,
+ bg = rep(c(0, 2), c(250, 10)))
+for ( a in a184 ) {
+ points(parameters(a), col = gray(0.7))
+ lines(ahull(a), col = gray(0.7))
+}
+
+# => For a given k, breakdown point 0. This means we can make the
+# corresponding archetype arbitrarily "large" just by changing
+# any of the data points.
+
+
+
+### One region with outliers: ########################################
+
+toy.outlier <- function(n, mean, sigma, ...) {
+ require(mvtnorm)
+ data(toy, package = 'archetypes')
+
+ for ( i in seq(length = length(n)) )
+ toy <- rbind(toy, rmvnorm(n[i], mean[[i]], sigma[[i]]))
+
+ toy
+}
+
+
+set.seed(1234)
+toy.o1 <- toy.outlier(5, mean = list(c(0, 30)),
+ sigma = list(diag(0.5, 2)))
+
+plot(toy.o1)
+
+
+
+### Original archetypes:
+
+set.seed(1234)
+a1 <- archetypes(toy.o1, 3)
+
+plot(a1, toy.o1)
+plot(a1, toy.o1, adata.show = TRUE)
+
+residuals.diagplot(a1)
+rss.diagplot(a1)
+
+
+
+### Weighted archetypes:
+
+w <- rep(c(1, 0), c(250, 5))
+
+set.seed(1234)
+wa1 <- archetypes(toy.o1, 3, family = archetypesFamily('weighted'), weights = w)
+
+plot(wa1, toy.o1)
+plot(wa1, toy.o1, adata.show = TRUE)
+
+
+
+### Robust archetypes:
+
+set.seed(1234)
+ra1 <- archetypes(toy.o1, 3, family = archetypesFamily('robust'))
+
+plot(ra1, toy.o1, adata.show = TRUE)
+
+residuals.diagplot(ra1)
+rss.diagplot(ra1)
+weights.diagplot(ra1, weights.type = 'reweights')
+
+par(mfrow = c(6, 7), mar = c(0, 0, 0, 0))
+movieplot(ra1, toy.o1, adata.show = TRUE, link.col.show = FALSE,
+ link.lty = 0, axes = FALSE, postfn = function(iter) box())
+
+# => looks like it does the right thing, but ...
+
+set.seed(1235)
+ra2 <- archetypes(toy.o1, 3, family = archetypesFamily('robust'))
+
+plot(ra2, toy.o1, adata.show = TRUE)
+
+weights.diagplot(ra2, weights.type = 'reweights')
+
+par(mfrow = c(1, 6), mar = c(0, 0, 0, 0))
+movieplot(ra2, toy.o1, adata.show = TRUE, link.col.show = FALSE,
+ link.lty = 0, axes = FALSE, postfn = function(iter) box())
+
+# => hmm ... it really depends on the starting values. What about
+# the weights?
+
+reweights.diagplot(ra1)
+reweights.diagplot(ra2)
+
+reweights.diagplot(ra1, highlight = c(1, 2))
+i.reweights.diagplot(ra1, c(1, 2))
+
+
+
+### Robust archetypes on a data set without outliers:
+
+set.seed(1234)
+a2 <- archetypes(toy, 3, family = archetypesFamily('robust'))
+
+plot(a2, toy, adata.show = TRUE)
+
+par(mfrow = c(6, 6), mar = c(0, 0, 0, 0))
+movieplot(a2, toy, adata.show = TRUE, link.col.show = FALSE,
+ link.lty = 0, axes = FALSE, postfn = function(iter) box())
+
+reweights.diagplot(a2)
+i.reweights.diagplot(a2, 1:50)
+
+
+
+# => Hmm ...
+
+set.seed(1234)
+a3 <- stepArchetypes(toy, family = archetypesFamily('robust'), k = 3, nrep = 5)
+
+plot(a3, toy)
+
+# => Not the expected behavior!
+
+
Property changes on: branches/pkg-robust/sandbox/demo.R
___________________________________________________________________
Name: svn:keywords
+ Date Revision Author URL Id
Name: svn:eol-style
+ native
More information about the Archetypes-commits
mailing list