[Archetypes-commits] r26 - in branches/pkg-robust: . R sandbox
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Feb 23 17:30:42 CET 2010
Author: manuel
Date: 2010-02-23 17:30:42 +0100 (Tue, 23 Feb 2010)
New Revision: 26
Modified:
branches/pkg-robust/DESCRIPTION
branches/pkg-robust/R/archetypes-diagplots.R
branches/pkg-robust/R/archetypes-kit-blocks.R
branches/pkg-robust/sandbox/demo.R
Log:
Modified: branches/pkg-robust/DESCRIPTION
===================================================================
--- branches/pkg-robust/DESCRIPTION 2010-02-22 17:20:18 UTC (rev 25)
+++ branches/pkg-robust/DESCRIPTION 2010-02-23 16:30:42 UTC (rev 26)
@@ -3,7 +3,7 @@
Title: Archetypal Analysis
Version: 1.0
Date: 2009-04-23
-Depends: nnls (>= 1.1)
+Depends: nnls (>= 1.1), modeltools
Suggests: MASS, vcd
Author: Manuel J. A. Eugster <manuel.eugster at stat.uni-muenchen.de>
Maintainer: Manuel J. A. Eugster <manuel.eugster at stat.uni-muenchen.de>
Modified: branches/pkg-robust/R/archetypes-diagplots.R
===================================================================
--- branches/pkg-robust/R/archetypes-diagplots.R 2010-02-22 17:20:18 UTC (rev 25)
+++ branches/pkg-robust/R/archetypes-diagplots.R 2010-02-23 16:30:42 UTC (rev 26)
@@ -72,32 +72,59 @@
-reweights.diagplot <- function(object, highlight = NULL,
- highlight.col = (seq(length(highlight)) + 1), ...) {
+reweights.diagplot <- function(object, col = 1, pch = 1, highlight = NULL,
+ highlight.col = (seq(length(highlight)) + 1),
+ highlight.pch = 13, ...) {
y <- rev(lapply(object$history, function(x) x[[1]]$reweights))
x <- seq(along = y[[1]])
- col <- rep(1, length = length(x))
+
+ col <- rep(col, length = length(x))
col[highlight] <- highlight.col
+ pch <- rep(pch, length = length(x))
+ pch[highlight] <- highlight.pch
+
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', ...)
+ plot(x, y[[i]], type = 'p', col = col, pch = pch, ylim = c(0, 1),
+ xlab = 'Index', ylab = 'Reweights', ...)
}
-i.reweights.diagplot <- function(object, i, lty = 1,
- col = (seq(length(i)) + 1), ...) {
+reweights.curve.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', ...)
+ xlab = 'Iterations', ylim = c(0, 1), ...)
}
+reweights.liftoff.diagplot <- function(object, ...) {
+
+ y <- sapply(object$history, function(x) x[[1]]$reweights)
+ y <- rev(colSums(y != 0))
+
+ barplot(y, xlab = 'Iterations', ylab = 'Reweights > 0', ...)
+}
+
+
+
+reweights.rss.diagplot <- function(object, ...) {
+ y1 <- sapply(object$history, function(x) x[[1]]$reweights)
+ y1 <- rev(colSums(y1))
+
+ y2 <- rev(sapply(object$history, function(x) x[[1]]$rss))
+
+ x <- seq(along = y1)
+
+ par(mfrow = c(2, 1))
+ plot(x, y1, type = 'l', xlab = 'Iterations', ylab = 'Reweights', ...)
+ plot(x, y2, type = 'l', xlab = 'Iterations', ylab = 'RSS', ...)
+}
Modified: branches/pkg-robust/R/archetypes-kit-blocks.R
===================================================================
--- branches/pkg-robust/R/archetypes-kit-blocks.R 2010-02-22 17:20:18 UTC (rev 25)
+++ branches/pkg-robust/R/archetypes-kit-blocks.R 2010-02-23 16:30:42 UTC (rev 26)
@@ -297,13 +297,8 @@
### Reweighting functions:
-#'
-#'
bisquare0.reweightsfn <- function(resid) {
- if ( is.null(resid) )
- return(NULL)
-
- resid <- apply(resid, 2, function(.) sum(abs(.)))
+ resid <- apply(resid, 2, function(x) sum(abs(x)))
resid0 <- resid < sqrt(.Machine$double.eps)
s <- resid / 6 * median(resid[!resid0])
@@ -311,11 +306,19 @@
ifelse(s < 1, (1 - s^2)^2, 0)
}
+bisquare.reweightsfn <- function(resid) {
+ resid.abs <- apply(resid, 2, function(x) sum(abs(x)))
+
+ mar <- mad(resid.abs, constant = 1) / 0.6754
+ k <- 4.685 * mar
+
+ resid.euc <- apply(resid / k, 2, function(x) sum(x^2))
+
+ ifelse(resid.abs <= k, (1 - resid.euc)^2, 0)
+}
+
tricube.reweightsfn <- function(resid) {
- if ( is.null(resid) )
- return(NULL)
-
- resid <- apply(resid, 2, function(.) sum(abs(.)))
+ resid <- apply(resid, 2, function(x) sum(abs(x)))
ifelse(resid < 1, (1 - resid^3)^3, 0)
}
Modified: branches/pkg-robust/sandbox/demo.R
===================================================================
--- branches/pkg-robust/sandbox/demo.R 2010-02-22 17:20:18 UTC (rev 25)
+++ branches/pkg-robust/sandbox/demo.R 2010-02-23 16:30:42 UTC (rev 26)
@@ -1,8 +1,14 @@
library(archetypes.dev.robust)
+library(modeltools)
+sapply(list.files('../R', full = TRUE), source, echo = TRUE)
+load('../data/toy.RData')
+
+
+
### Archetypes with one outlier; breakdown point simulations ########
data(toy)
@@ -112,6 +118,7 @@
plot(ra2, toy.o1, adata.show = TRUE)
weights.diagplot(ra2, weights.type = 'reweights')
+rss.diagplot(ra2)
par(mfrow = c(1, 6), mar = c(0, 0, 0, 0))
movieplot(ra2, toy.o1, adata.show = TRUE, link.col.show = FALSE,
@@ -124,10 +131,16 @@
reweights.diagplot(ra2)
reweights.diagplot(ra1, highlight = c(1, 2))
-i.reweights.diagplot(ra1, c(1, 2))
+reweights.curve.diagplot(ra1, c(1, 2))
+reweights.liftoff.diagplot(ra1)
+reweights.liftoff.diagplot(ra2)
+reweights.rss.diagplot(ra1)
+reweights.rss.diagplot(ra2)
+
+
### Robust archetypes on a data set without outliers:
set.seed(1234)
@@ -140,10 +153,8 @@
link.lty = 0, axes = FALSE, postfn = function(iter) box())
reweights.diagplot(a2)
-i.reweights.diagplot(a2, 1:50)
+reweights.curve.diagplot(a2, 1:50)
-
-
# => Hmm ...
set.seed(1234)
@@ -151,6 +162,23 @@
plot(a3, toy)
-# => Not the expected behavior!
+# => Not the intended behavior!
+
+### So, as always, it is not the easy solution ... ###################
+
+set.seed(1235)
+ra2 <- archetypes(toy.o1, 3, family = archetypesFamily('robust',
+ reweightsfn = bisquare.reweightsfn))
+
+plot(ra2, toy.o1, adata.show = TRUE)
+
+weights.diagplot(ra2, weights.type = 'reweights')
+rss.diagplot(ra2)
+
+par(mfrow = c(2, 5), 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())
+
+reweights.rss.diagplot(ra2)
More information about the Archetypes-commits
mailing list