[Archetypes-commits] r32 - in branches/pkg-robust: R sandbox
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Mar 18 15:53:06 CET 2010
Author: manuel
Date: 2010-03-18 15:53:06 +0100 (Thu, 18 Mar 2010)
New Revision: 32
Added:
branches/pkg-robust/R/archetypes-panorama.R
branches/pkg-robust/R/archetypes-robust.R
branches/pkg-robust/R/archetypes-weighted.R
branches/pkg-robust/sandbox/Ozone.R
branches/pkg-robust/sandbox/demo2.R
Modified:
branches/pkg-robust/R/archetypes-barplot.R
branches/pkg-robust/R/archetypes-class.R
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
branches/pkg-robust/sandbox/demo.R
Log:
Modified: branches/pkg-robust/R/archetypes-barplot.R
===================================================================
--- branches/pkg-robust/R/archetypes-barplot.R 2010-03-12 10:19:34 UTC (rev 31)
+++ branches/pkg-robust/R/archetypes-barplot.R 2010-03-18 14:53:06 UTC (rev 32)
@@ -3,7 +3,10 @@
#' Barplot of archetypes.
#' @param height An \code{\link{archetypes}} object.
#' @param data The original data matrix.
-#' @param beside Portray the archetypes as stacked bars.
+#' @param which \code{below} creates a barplot for each archetype,
+#' \code{beside} creates one barplot with bars side by side.
+#' @param which.beside Barplot according to \code{atypes} or \code{variables}.
+#' @param which.below \code{compressed} plots the labels only once.
#' @param percentage Show real values or percentages according to the
#' original data.
#' @param ... Passed to the underlying \code{\link{barplot}} call.
@@ -11,30 +14,86 @@
#' @method barplot archetypes
#' @importFrom graphics barplot
#' @export
-barplot.archetypes <- function(height, data, beside=TRUE, percentage=FALSE, ...) {
+barplot.archetypes <- function(height, data,
+ which = c('below', 'beside'),
+ which.beside = c('atypes', 'variables'),
+ which.below = c('compressed', 'default'),
+ percentiles = FALSE, ...) {
+
+ ### Helpers:
+ .beside.atypes <- function() {
+ barplot(t(atypes), ylab=ylab, beside=TRUE, ylim=ylim, ...)
+ }
+
+
+ .beside.variables <- function() {
+ barplot(atypes, ylab=ylab, beside=TRUE, ylim=ylim, ...)
+ }
+
+
+ .below.default <- function() {
+ p <- nrow(atypes)
+
+ layout(matrix(1:p, nrow = p, byrow = TRUE))
+ for ( i in 1:p )
+ barplot(atypes[i,], main=paste('Archetype', i),
+ ylab=ylab, ylim=ylim, ...)
+ }
+
+
+ .below.compressed <- function() {
+ p <- nrow(atypes) + 1
+
+ layout(matrix(1:p, nrow = p, byrow = TRUE))
+ for ( i in 1:(p - 1) ) {
+ par(mar = c(0, 5, 1, 0) + 0.1)
+ x.at <- barplot(atypes[i,], ylab=ylab, ylim=ylim,
+ names.arg='', las=2, ...)
+ mtext(sprintf('Archetype %s', i), side = 2, line = 4,
+ cex = par('cex'))
+ }
+ text(x.at, par("usr")[3] - 1, srt = 90, adj = 1,
+ labels = colnames(atypes), xpd = NA)
+ }
+
+
+ .perc <- function(x, data, digits = 0) {
+ Fn <- ecdf(data)
+ round(Fn(x) * 100, digits = digits)
+ }
+
+
+ ### Plot:
+ opar <- par(no.readonly = TRUE)
+ on.exit(par(opar))
+
+ which <- match.arg(which)
+
+ if ( which == 'beside' )
+ which.arg <- match.arg(which.beside)
+ else
+ which.arg <- match.arg(which.below)
+
atypes <- atypes(height)
+ rownames(atypes) <- sprintf('Archetype %s',
+ seq(length = nrow(atypes)))
- if ( !percentage ) {
+ if ( !percentiles ) {
ylab <- 'Value'
ylim <- NULL
}
else {
- m <- sapply(data, max)
- atypes <- t(t(atypes) / m * 100)
- ylab <- 'Percentage'
- ylim <- c(0,100)
- }
+ atypes <- sapply(seq(length = ncol(data)),
+ function(i)
+ .perc(atypes[, i], data[, i]))
+ colnames(atypes) <- colnames(data)
- if ( beside ) {
- barplot(atypes, ylab=ylab, beside=TRUE, ylim=ylim, ...)
+ ylab <- 'Percentile'
+ ylim <- c(0, 100)
}
- else {
- p <- nrow(atypes)
-
- par(mfrow=c(p,1))
- for ( i in 1:p )
- barplot(atypes[i,], main=paste('Archetype', i),
- ylab=ylab, ylim=ylim, ...)
- }
+
+ do.call(sprintf('.%s.%s', which, which.arg), list())
+
+ invisible(atypes)
}
Modified: branches/pkg-robust/R/archetypes-class.R
===================================================================
--- branches/pkg-robust/R/archetypes-class.R 2010-03-12 10:19:34 UTC (rev 31)
+++ branches/pkg-robust/R/archetypes-class.R 2010-03-18 14:53:06 UTC (rev 32)
@@ -15,8 +15,10 @@
#' @param betas The data coefficients; a \eqn{p \times n} matrix.
#' @param zas The temporary archetypes.
#' @param family The archetypes family.
+#' @param familyArgs Additional arguments for family blocks.
#' @param residuals The residuals.
#' @param weights The data weights.
+#' @param reweights The data reweights.
#' @return A list with an element for each parameter and class attribute
#' \code{archetypes}.
#' @seealso \code{\link{archetypes}}, \code{\link{atypes}}, \code{\link{ntypes}},
@@ -25,8 +27,8 @@
#' @export
as.archetypes <- function(archetypes, k, alphas, rss, iters = NULL, call = NULL,
history = NULL, kappas = NULL, betas = NULL, zas = NULL,
- family = NULL, residuals = NULL, weights = NULL,
- reweights = NULL) {
+ family = NULL, familyArgs = NULL, residuals = NULL,
+ weights = NULL, reweights = NULL) {
return(structure(list(archetypes = archetypes,
k = k,
@@ -39,6 +41,7 @@
call = call,
history = history,
family = family,
+ familyArgs = familyArgs,
residuals = residuals,
weights = weights,
reweights = reweights),
@@ -46,8 +49,6 @@
}
setOldClass('archetypes')
-setOldClass('weightedArchetypes')
-setOldClass('robustArchetypes')
@@ -94,8 +95,6 @@
#' @importFrom modeltools parameters
setMethod('parameters', 'archetypes', parameters.archetypes)
-setMethod('parameters', 'weightedArchetypes', parameters.archetypes)
-setMethod('parameters', 'robustArchetypes', parameters.archetypes)
Modified: branches/pkg-robust/R/archetypes-diagplots.R
===================================================================
--- branches/pkg-robust/R/archetypes-diagplots.R 2010-03-12 10:19:34 UTC (rev 31)
+++ branches/pkg-robust/R/archetypes-diagplots.R 2010-03-18 14:53:06 UTC (rev 32)
@@ -1,4 +1,3 @@
-### distanz von jedem archetyp zum nächsten punkt
residuals.diagplot <- function(object, ref.order = 1, ...) {
y <- residuals(object)
@@ -143,14 +142,54 @@
if ( is.null(ref.order) )
ix <- x
else
- ix <- order(d[, 1])
+ ix <- order(d[, ref.order])
ylim <- c(0, max(d))
par(mfrow = c(ncol(d), 1))
for ( i in seq(length = ncol(d)) )
- plot(x, d[ix, i], xlab = sprintf('Archetype %s', i), ylab = 'Distance',
- ylim = ylim, ...)
+ plot(x, d[ix, i], xlab = sprintf('Archetype %s', i),
+ ylab = 'Distance', ylim = ylim, ...)
invisible(d)
}
+
+
+archetypes.distance.diagplot <- function(object, data,
+ distfn = distEuclidean, ...) {
+
+ opar <- par()
+
+ d <- distfn(data, parameters(object))
+ x <- seq(length = nrow(d))
+
+ ylim <- c(0, max(d))
+
+ globals <- array(dim = c(dim(d), 2))
+
+ par(mfrow = c(ncol(d) + 1, 1), xpd = NA)
+ for ( i in seq(length = ncol(d)) ) {
+ ix <- order(d[, i])
+
+ plot(x, d[ix, i], xlab = '', ylab = 'Distance', ylim = ylim,
+ type = 'b', axes = FALSE, ...)
+ axis(2)
+ box()
+
+ globals[, i, 1] <- grconvertX(rank(d[, i]), to = 'device')
+ globals[, i, 2] <- grconvertY(d[, i], to = 'device')
+ }
+ axis(1)
+ mtext('Data index', side = 1, line = 3, cex = par('cex'))
+
+ for ( i in seq(length = nrow(globals)) ) {
+ px <- grconvertX(globals[i, , 1], from = 'device')
+ py <- grconvertY(globals[i, , 2], from = 'device')
+
+ lines(px, py, col = gray(0.7))
+ }
+
+ par(opar)
+
+ invisible(d)
+}
Modified: branches/pkg-robust/R/archetypes-kit-blocks.R
===================================================================
--- branches/pkg-robust/R/archetypes-kit-blocks.R 2010-03-12 10:19:34 UTC (rev 31)
+++ branches/pkg-robust/R/archetypes-kit-blocks.R 2010-03-18 14:53:06 UTC (rev 32)
@@ -5,7 +5,7 @@
#' Scaling block: standardize to mean 0 and standard deviation 1.
#' @param x Data matrix.
#' @return Standardized data matrix with some attribues.
-std.scalefn <- function(x) {
+std.scalefn <- function(x, ...) {
m = rowMeans(x)
x = x - m
@@ -21,7 +21,7 @@
#' @param x Standardized data matrix.
#' @param zs Archetypes matrix
#' @return Rescaled archetypes.
-std.rescalefn <- function(x, zs) {
+std.rescalefn <- function(x, zs, ...) {
m = attr(x, '.Meta')$mean
s = attr(x, '.Meta')$sd
@@ -37,7 +37,7 @@
#' Scaling block: no scaling.
#' @param x Data matrix.
#' @return Data matrix.
-no.scalefn <- function(x) {
+no.scalefn <- function(x, ...) {
return(x)
}
@@ -45,7 +45,7 @@
#' @param x Data matrix.
#' @param zs Archetypes matrix.
#' @return Archetypes zs.
-no.rescalefn <- function(x, zs) {
+no.rescalefn <- function(x, zs, ...) {
return(zs)
}
@@ -60,7 +60,7 @@
#' with an additonal row containing \code{huge} values.
make.dummyfn <- function(huge=200) {
- bp.dummyfn <- function(x) {
+ bp.dummyfn <- function(x, ...) {
y = rbind(x, rep(huge, ncol(x)))
attr(y, '.Meta') = attr(x, '.Meta')
@@ -77,7 +77,7 @@
#' @param x Data matrix.
#' @param zs Archetypes matrix.
#' @return Archetypes zs.
-rm.undummyfn <- function(x, zs) {
+rm.undummyfn <- function(x, zs, ...) {
dr = attr(x, '.Meta')$dummyrow
return(zs[-dr,])
@@ -87,7 +87,7 @@
#' Dummy block: no dummy row.
#' @param x Data matrix.
#' @return Data matrix x.
-no.dummyfn <- function(x) {
+no.dummyfn <- function(x, ...) {
return(x)
}
@@ -95,7 +95,7 @@
#' @param x Data matrix.
#' @param zs Archetypes matrix.
#' @return Archetypes zs.
-no.undummyfn <- function(x, zs) {
+no.undummyfn <- function(x, zs, ...) {
return(zs)
}
@@ -107,7 +107,7 @@
#' @param alphas The coefficients.
#' @param x Data matrix.
#' @return The solved linear system.
-qrsolve.zalphasfn <- function(alphas, x) {
+qrsolve.zalphasfn <- function(alphas, x, ...) {
return(t(qr.solve(alphas %*% t(alphas)) %*% alphas %*% t(x)))
}
@@ -117,7 +117,7 @@
#' @param alphas The coefficients.
#' @param x Data matrix.
#' @return The solved linear system.
-ginv.zalphasfn <- function(alphas, x) {
+ginv.zalphasfn <- function(alphas, x, ...) {
require(MASS)
return(t(ginv(alphas %*% t(alphas)) %*% alphas %*% t(x)))
@@ -129,7 +129,7 @@
#' @param alphas The coefficients.
#' @param x Data matrix.
#' @return The solved linear system.
-opt.zalphasfn <- function(alphas, x) {
+opt.zalphasfn <- function(alphas, x, ...) {
z <- rnorm(nrow(x)*nrow(alphas))
fun <- function(z){
@@ -152,7 +152,7 @@
#' @param C The archetypes matrix.
#' @param d The data matrix.
#' @return Recalculated alpha.
-nnls.alphasfn <- function(coefs, C, d) {
+nnls.alphasfn <- function(coefs, C, d, ...) {
require(nnls)
n = ncol(d)
@@ -168,7 +168,7 @@
#' @param C The archetypes matrix.
#' @param d The data matrix.
#' @return Recalculated alpha.
-snnls.alphasfn <- function(coefs, C, d) {
+snnls.alphasfn <- function(coefs, C, d, ...) {
require(nnls)
n = ncol(d)
@@ -215,7 +215,7 @@
#' Norm block: standard matrix norm (spectral norm).
#' @param m Matrix.
#' @return The norm.
-norm2.normfn <- function(m) {
+norm2.normfn <- function(m, ...) {
return(max(svd(m)$d))
}
@@ -223,7 +223,7 @@
#' Norm block: euclidian norm.
#' @param m Matrix.
#' @return The norm.
-euc.normfn <- function(m) {
+euc.normfn <- function(m, ...) {
return(sum(apply(m, 2, function(x){sqrt(sum(x^2))})))
}
@@ -236,7 +236,7 @@
#' @return A function which returns a list with alpha and beta.
make.random.initfn <- function(k) {
- bp.initfn <- function(x, p) {
+ bp.initfn <- function(x, p, ...) {
n = ncol(x)
b = matrix(0, nrow=n, ncol=p)
@@ -257,7 +257,7 @@
#' @return A function which returns a list with alpha and beta.
make.fix.initfn <- function(indizes) {
- fix.initfn <- function(x, p) {
+ fix.initfn <- function(x, p, ...) {
n = ncol(x)
b = matrix(0, nrow = n, ncol = p)
@@ -279,17 +279,23 @@
#' @param data A numeric \eqn{m \times n} data matrix.
#' @param weights Vector of data weights within \eqn{[0, 1]}.
#' @return Weighted data matrix.
-center.weightfn <- function(data, weights) {
+center.weightfn <- function(data, weights, ...) {
if ( is.null(weights) )
return(data)
+ weights <- as.numeric(1 - weights)
+
dr <- attr(data, '.Meta')$dummyrow
- weights <- as.numeric(1 - weights)
- center <- rowMeans(data[-dr,])
+ if ( is.null(dr) ) {
+ center <- rowMeans(data)
+ data <- data + t(weights * t(center - data))
+ }
+ else {
+ center <- rowMeans(data[-dr, ])
+ data[-dr, ] <- data[-dr, ] + t(weights * t(center - data[-dr, ]))
+ }
- data[-dr,] <- data[-dr,] + t(weights * t(center - data[-dr,]))
-
data
}
@@ -297,16 +303,25 @@
### Reweighting functions:
-bisquare0.reweightsfn <- function(resid, reweights) {
+bisquare0.reweightsfn <- function(resid, reweights, ...) {
resid <- apply(resid, 2, function(x) sum(abs(x)))
resid0 <- resid < sqrt(.Machine$double.eps)
- s <- resid / 6 * median(resid[!resid0])
+ s <- 6 * median(resid[!resid0])
+ v <- resid / s
- ifelse(s < 1, (1 - s^2)^2, 0)
+ ifelse(v < 1, (1 - v^2)^2, 0)
}
-bisquare.reweightsfn <- function(resid, reweights) {
+
+binary.bisquare0.reweightsfn <- function(resid, reweights,
+ cutpoint = 0.1, ...) {
+ rw <- bisquare0.reweightsfn(resid, reweights, ...)
+ ifelse(rw < cutpoint, 0, 1)
+}
+
+
+bisquare.reweightsfn <- function(resid, reweights, ...) {
resid.abs <- apply(resid, 2, function(x) sum(abs(x)))
mar <- mad(resid.abs, constant = 1) / 0.6754
@@ -317,13 +332,37 @@
ifelse(resid.abs <= k, (1 - resid.euc)^2, 0)
}
-tricube.reweightsfn <- function(resid, reweights) {
+tricube.reweightsfn <- function(resid, reweights, ...) {
resid <- apply(resid, 2, function(x) sum(abs(x)))
ifelse(resid < 1, (1 - resid^3)^3, 0)
}
+leastwise.tricube.reweightsfn <- function(resid, reweights,
+ nperc = 0.8, minrw = 0.3, ...) {
+ resid <- apply(resid, 2, function(x) sum(abs(x)))
+ resid0 <- resid < sqrt(.Machine$double.eps)
+ resido <- order(resid)
+ n0 <- ceiling(sum(!resid0) * nperc)
+
+ ix <- resido[!resid0][seq(length = n0)]
+
+ rw <- ifelse(resid < 1, (1 - resid^3)^3, 0)
+ rw[ix] <- ifelse(rw[ix] < minrw, minrw, rw[ix])
+
+ rw
+}
+
+
+binary.tricube.reweightsfn <- function(resid, reweights,
+ cutpoint = 0.1, ...) {
+ rw <- tricube.reweightsfn(resid, reweights, ...)
+ ifelse(rw < cutpoint, 0, 1)
+}
+
+
+
### Archetypes family:
#' Archetypes family constructor.
@@ -373,18 +412,3 @@
class = NULL)
}
-.weighted.archetypesFamily <- function() {
- f <- .original.archetypesFamily()
- f$class <- 'weightedArchetypes'
- f$weightfn <- center.weightfn
- f
-}
-
-.robust.archetypesFamily <- function() {
- f <- .original.archetypesFamily()
- f$class <- 'robustArchetypes'
- f$weightfn <- center.weightfn
- f$reweightsfn <- tricube.reweightsfn
- f
-}
-
Modified: branches/pkg-robust/R/archetypes-kit.R
===================================================================
--- branches/pkg-robust/R/archetypes-kit.R 2010-03-12 10:19:34 UTC (rev 31)
+++ branches/pkg-robust/R/archetypes-kit.R 2010-03-18 14:53:06 UTC (rev 32)
@@ -7,6 +7,7 @@
#' Perform archetypal analysis on a data matrix.
#' @param data A numeric \eqn{n \times m} data matrix.
#' @param k The number of archetypes.
+#' @param weights Data weights vector.
#' @param maxIterations The maximum number of iterations.
#' @param minImprovement The minimal value of improvement between two
#' iterations.
@@ -16,6 +17,7 @@
#' further analyses.
#' @param family Blocks defining the underlying problem solving mechanisms;
#' see \code{\link{archetypesFamily}}.
+#' @param ... Additional arguments for family blocks.
#' @return An object of class \code{\link{archetypes}}, see
#' \code{\link{as.archetypes}}.
#' @seealso \code{\link{stepArchetypes}}
@@ -26,13 +28,14 @@
#' a <- archetypes(toy, 3)
#' @export
#' @note Please see the vignette for a detailed explanation!
-archetypes <- function(data, k, maxIterations = 100,
+archetypes <- function(data, k, weights = NULL, maxIterations = 100,
minImprovement = sqrt(.Machine$double.eps),
maxKappa = 1000, verbose = TRUE, saveHistory = TRUE,
- family = archetypesFamily('original'), weights = NULL) {
+ family = archetypesFamily('original'), ...) {
### Helpers:
mycall <- match.call()
+ famargs <- list(...)
history <- NULL
snapshot <- function(i) {
@@ -53,9 +56,9 @@
### Data preparation:
x1 <- t(data)
- x1 <- family$scalefn(x1)
- x1 <- family$dummyfn(x1)
- x0 <- family$weightfn(x1, weights)
+ x1 <- family$scalefn(x1, ...)
+ x1 <- family$dummyfn(x1, ...)
+ x0 <- family$weightfn(x1, weights, ...)
x <- x0
n <- ncol(x)
@@ -63,7 +66,7 @@
### Initialization:
- init <- family$initfn(x, k)
+ init <- family$initfn(x, k, ...)
betas <- init$betas
alphas <- init$alphas
@@ -72,7 +75,7 @@
zs <- x %*% betas
resid <- zs %*% alphas - x
- rss <- family$normfn(resid) / n
+ rss <- family$normfn(resid, ...) / n
reweights <- rep(1, n)
@@ -94,30 +97,30 @@
tryCatch(while ( (i <= maxIterations) & (imp >= minImprovement) ) {
## Reweight data:
- reweights <- family$reweightsfn(resid, reweights)
- x <- family$weightfn(x0, reweights)
+ reweights <- family$reweightsfn(resid, reweights, ...)
+ x <- family$weightfn(x0, reweights, ...)
## Alpha's:
- alphas <- family$alphasfn(alphas, zs, x)
- zas <- family$zalphasfn(alphas, x)
- rss1 <- family$normfn(zas %*% alphas - x) / n
+ alphas <- family$alphasfn(alphas, zs, x, ...)
+ zas <- family$zalphasfn(alphas, x, ...)
+ rss1 <- family$normfn(zas %*% alphas - x, ...) / n
kappas[c('alphas', 'zas')] <- c(kappa(alphas), kappa(zas))
## Beta's:
- betas <- family$betasfn(betas, x, zas)
+ betas <- family$betasfn(betas, x, zas, ...)
zs <- x %*% betas
kappas[c('betas', 'zs')] <- c(kappa(betas), kappa(zs))
## Residuals, RSS and improvement:
- alphas0 <- family$alphasfn(alphas, zs, x0)
+ alphas0 <- family$alphasfn(alphas, zs, x0, ...)
resid <- zs %*% alphas0 - x0
- rss2 <- family$normfn(resid) / n
+ rss2 <- family$normfn(resid, ...) / n
imp <- rss - rss2
rss <- rss2
@@ -136,19 +139,20 @@
i <- i + 1
},
- error=function(e) errormsg <<- e)
+ error = function(e) errormsg <<- e)
### Check illness:
if ( !is.null(errormsg) ) {
warning('k=', k, ': ', errormsg)
- return(as.archetypes(NULL, k, NULL, NA, iters=i,
- call=mycall, history=history, kappas=kappas))
+ return(as.archetypes(NULL, k, NULL, NA, iters = i,
+ call = mycall, history = history,
+ kappas = kappas))
}
if ( any(isIll) )
- warning('k=', k, ': ', paste(names(isIll)[isIll], collapse=', '),
- ' > maxKappa', sep='')
+ warning('k=', k, ': ', paste(names(isIll)[isIll], collapse = ', '),
+ ' > maxKappa', sep = '')
### Rescale archetypes, etc.:
@@ -168,8 +172,8 @@
return(as.archetypes(t(zs), k, t(alphas), rss, iters = (i-1),
call = mycall, history = history, kappas = kappas,
betas = t(betas), family = family,
- residuals = t(resid), weights = weights,
- reweights = reweights))
+ familyArgs = famargs, residuals = t(resid),
+ weights = weights, reweights = reweights))
}
Modified: branches/pkg-robust/R/archetypes-movie.R
===================================================================
--- branches/pkg-robust/R/archetypes-movie.R 2010-03-12 10:19:34 UTC (rev 31)
+++ branches/pkg-robust/R/archetypes-movie.R 2010-03-18 14:53:06 UTC (rev 32)
@@ -15,22 +15,39 @@
#' ssleep=0, bsleep=0, postfn=function(iter){}, ...)
#' @aliases movieplot
#' @export
-movieplot <- function(zs, data, show=c('atypes', 'adata'),
- ssleep=0, bsleep=0, postfn=function(iter){}, ...) {
+movieplot <- function(zs, data, show=c('atypes', 'adata', 'rwdata'),
+ ssleep=0, bsleep=0, postfn=function(iter){},
+ rwdata.col1 = gray(0.7), rwdata.col2 = 2, ...) {
+ show <- match.arg(show)
steps <- length(zs$history)
- atypesmovie <- ifelse(show[1] == 'atypes', TRUE, FALSE)
+ if ( show == 'rwdata' )
+ data <- zs$family$scalefn(t(data))
+
+ # ... and play:
Sys.sleep(ssleep)
- # ... and play:
for ( i in seq_len(steps)-1 ) {
a <- ahistory(zs, step=i)
- if ( atypesmovie )
- plot(a, data, ...)
- else
- plot(adata(a), ...)
+ switch(show,
+
+ atypes = {
+ plot(a, data, ...)
+ },
+
+ adata = {
+ plot(adata(a), ...)
+ },
+
+ rwdata = {
+ d <- zs$family$weightfn(data, a$reweights)
+
+ plot(t(data), col = rwdata.col1, ...)
+ points(t(d), col = rwdata.col2, ...)
+ })
+
postfn(i)
Sys.sleep(bsleep)
Added: branches/pkg-robust/R/archetypes-panorama.R
===================================================================
--- branches/pkg-robust/R/archetypes-panorama.R (rev 0)
+++ branches/pkg-robust/R/archetypes-panorama.R 2010-03-18 14:53:06 UTC (rev 32)
@@ -0,0 +1,55 @@
+
+
+panorama <- function(object, ...) {
+ UseMethod('panorama')
+}
+
+
+
+panorama.archetypes <- function(object, data, distfn = distEuclidean,
+ ref.order = NULL, xlab = 'Index',
+ ylab = 'Distance', col = 1, pch = 1,
+ atypes.col = (seq(length = ntypes(object)) + 1),
+ atypes.pch = rep(19, ntypes(object)), ...) {
+
+ n1 <- nrow(data)
+ n2 <- ntypes(object)
+
+ data <- rbind(data, parameters(object))
+ dist <- distfn(data, parameters(object))
+
+ x <- seq(length = n1 + n2)
+
+ col <- c(rep(col, n1), atypes.col)
+ pch <- c(rep(pch, n1), atypes.pch)
+
+ if ( is.null(ref.order) )
+ ix <- x
+ else
+ ix <- order(dist[, ref.order])
+
+
+ opar <- par(no.readonly = TRUE)
+ on.exit(par(opar))
+
+ ylim <- c(0, max(dist))
+
+ mar <- opar$mar
+ mar[2] <- max(mar[2], 5)
+
+ par(mfrow = c(n2, 1), mar = mar)
+
+ for ( i in seq(length = n2) ) {
+ plot(x, dist[ix, i], ylim = ylim, xlab = xlab,
+ ylab = ylab, col = col[ix], pch = pch[ix], ...)
+
+ mtext(sprintf('Archetype %s', i), side = 2, line = 4,
+ cex = par('cex'))
+ }
+
+
+ invisible(dist)
+}
+
+
+
Property changes on: branches/pkg-robust/R/archetypes-panorama.R
___________________________________________________________________
Name: svn:keywords
+ Date Revision Author URL Id
Name: svn:eol-style
+ native
Modified: branches/pkg-robust/R/archetypes-plot.R
===================================================================
--- branches/pkg-robust/R/archetypes-plot.R 2010-03-12 10:19:34 UTC (rev 31)
+++ branches/pkg-robust/R/archetypes-plot.R 2010-03-18 14:53:06 UTC (rev 32)
@@ -81,7 +81,6 @@
}
-
#' Plot weighted archetypes.
plot.weightedArchetypes <- function(x, y,
adata.show = FALSE, data.col = 1, data.pch = 21,
Added: branches/pkg-robust/R/archetypes-robust.R
===================================================================
--- branches/pkg-robust/R/archetypes-robust.R (rev 0)
+++ branches/pkg-robust/R/archetypes-robust.R 2010-03-18 14:53:06 UTC (rev 32)
@@ -0,0 +1,24 @@
+#' @include archetypes-kit-blocks.R
+{}
+
+
+
+### Class and Family:
+
+.robust.archetypesFamily <- function() {
+ f <- .original.archetypesFamily()
+ f$class <- 'robustArchetypes'
+ f$weightfn <- center.weightfn
+ f$reweightsfn <- bisquare0.reweightsfn
+ f
+}
+
+setOldClass('robustArchetypes')
+
+
+
+### Methods:
+
+#' @importFrom modeltools parameters
+setMethod('parameters', 'robustArchetypes', parameters.archetypes)
+
Property changes on: branches/pkg-robust/R/archetypes-robust.R
___________________________________________________________________
Name: svn:keywords
+ Date Revision Author URL Id
Name: svn:eol-style
+ native
Added: branches/pkg-robust/R/archetypes-weighted.R
===================================================================
--- branches/pkg-robust/R/archetypes-weighted.R (rev 0)
+++ branches/pkg-robust/R/archetypes-weighted.R 2010-03-18 14:53:06 UTC (rev 32)
@@ -0,0 +1,24 @@
+#' @include archetypes-kit-blocks.R
+{}
+
+
+
+### Class and Family:
+
+.weighted.archetypesFamily <- function() {
+ f <- .original.archetypesFamily()
+ f$class <- 'weightedArchetypes'
+ f$weightfn <- center.weightfn
+ f
+}
+
+setOldClass('weightedArchetypes')
+
+
+
+### Methods:
+
+#' @importFrom modeltools parameters
+setMethod('parameters', 'weightedArchetypes', parameters.archetypes)
+
+
Property changes on: branches/pkg-robust/R/archetypes-weighted.R
___________________________________________________________________
Name: svn:keywords
+ Date Revision Author URL Id
Name: svn:eol-style
+ native
Added: branches/pkg-robust/sandbox/Ozone.R
===================================================================
--- branches/pkg-robust/sandbox/Ozone.R (rev 0)
+++ branches/pkg-robust/sandbox/Ozone.R 2010-03-18 14:53:06 UTC (rev 32)
@@ -0,0 +1,87 @@
+
+#library(archetypes.dev.robust)
+
+library(flexclust)
+library(modeltools)
+sapply(list.files('../R', full = TRUE), source, echo = TRUE)
+
+
+
+### Ozone data:
+
+data('Ozone', package = 'mlbench')
+
+oz <- Ozone[, -c(1, 2, 3, 9)]
+oz <- na.omit(oz)
+colnames(oz) <- c('OZONE', '500MH', 'WDSP', 'HMDTY', 'STMP',
+ 'INVHT', 'PRGRT', 'INVTMP', 'VZBLTY')
+
+oz0 <- scale(oz)
+
+
+
+### The three original archetypes:
+
+set.seed(1234)
+a.oz <- archetypes(oz0, 3)
+
+barplot(a.oz, oz0, percentiles = TRUE)
+
+panorama(a.oz, oz0)
+panorama(a.oz, oz0, ref.order = 1)
+
+
+
+### The robust archetypes:
+
+set.seed(1236)
+ra.oz <- archetypes(oz0, 3, family = archetypesFamily('robust'))
+
+weights(ra.oz, type = 'reweights')
+
+barplot(ra.oz, oz0, percentiles = TRUE)
+
+
+
+### Data set with outliers:
+
+outliers <- t(sapply(runif(5, min = 1.5, max = 2),
+ function(x)
+ x * apply(oz, 2, max) + apply(oz, 2, IQR)))
+
+oz1 <- scale(rbind(oz, outliers))
+
+
+
+### Original archetypes:
+
+set.seed(1234)
+a.oz1 <- archetypes(oz1, 3)
+
+barplot(a.oz, oz1, percentiles = TRUE)
+
+panorama(a.oz1, oz1)
+
+
+
+### Robust archetypes:
+
+set.seed(1236)
+ra.oz1 <- archetypes(oz1, 3, family = archetypesFamily('robust'))
+
+weights(ra.oz1, type = 'reweights')
+
+barplot(ra.oz1, oz1, percentiles = TRUE)
+
+panorama(ra.oz1, oz1)
+
+
+
+## Compared to the original archetypes:
+
+parameters(a.oz)
+parameters(ra.oz1)
+
+barplot(ra.oz1, oz1, percentiles = TRUE)
+x11(); barplot(a.oz, oz0, percentiles = TRUE)
+
Property changes on: branches/pkg-robust/sandbox/Ozone.R
___________________________________________________________________
Name: svn:keywords
+ Date Revision Author URL Id
Name: svn:eol-style
+ native
Modified: branches/pkg-robust/sandbox/demo.R
===================================================================
--- branches/pkg-robust/sandbox/demo.R 2010-03-12 10:19:34 UTC (rev 31)
+++ branches/pkg-robust/sandbox/demo.R 2010-03-18 14:53:06 UTC (rev 32)
@@ -84,9 +84,12 @@
rss.diagplot(a1, sort = TRUE)
archetypes.view.diagplot(a1, toy.o1)
+archetypes.view.diagplot(a1, toy.o1, ref.order = 1)
+archetypes.distance.diagplot(a1, toy.o1)
+
### Weighted archetypes:
w <- rep(c(1, 0), c(250, 5))
@@ -113,7 +116,7 @@
archetypes.view.diagplot(ra1, toy.o1)
archetypes.view.diagplot(ra1, toy.o1, ref.order = 1)
-par(mfrow = c(6, 7), mar = c(0, 0, 0, 0))
+par(mfrow = c(7, 6), 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())
@@ -155,13 +158,15 @@
plot(a2, toy, adata.show = TRUE)
-par(mfrow = c(6, 6), mar = c(0, 0, 0, 0))
+par(mfrow = c(5, 4), 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)
+reweights.liftoff.diagplot(a2)
reweights.curve.diagplot(a2, 1:50)
+
# => Hmm ...
set.seed(1234)
Added: branches/pkg-robust/sandbox/demo2.R
===================================================================
--- branches/pkg-robust/sandbox/demo2.R (rev 0)
+++ branches/pkg-robust/sandbox/demo2.R 2010-03-18 14:53:06 UTC (rev 32)
@@ -0,0 +1,176 @@
+
+#library(archetypes.dev.robust)
+
+library(flexclust)
+library(modeltools)
+sapply(list.files('../R', full = TRUE), source, echo = TRUE)
+
+
+
+### Data:
+
+load('../data/toy.RData')
+
+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)
+
+
+
+### Robust archetypes on a data set with outliers:
+
+set.seed(1234)
+ra1 <- archetypes(toy.o1, 3, family = archetypesFamily('robust'))
+
+plot(ra1, toy.o1, adata.show = TRUE)
+
+weights(ra1, type = 'reweights')
+
+
+## Stability:
+
+set.seed(1235)
+sra1 <- stepArchetypes(toy.o1, family = archetypesFamily('robust', reweightsfn = bisquare0.reweightsfn),
+ k = 3, nrep = 5)
+
+plot(sra1, toy.o1)
+
+
+
+### Robust archetypes on a data set without outliers:
+
+set.seed(1235)
+a2 <- archetypes(toy, 3, family = archetypesFamily('robust'))
+
+plot(a2, toy, adata.show = TRUE)
+
+weights(a2, type = 'reweights')
+
+
+## Stability:
+
+set.seed(1235)
+sa1 <- stepArchetypes(toy, family = archetypesFamily('robust'),
+ k = 3, nrep = 5)
+
+plot(sa1, toy)
+
+
+
+### Bisquare0 reweights function:
+
+## ... data set with outliers:
+
+opar <- par(mfrow = c(5, 3), mar = c(0, 0, 0, 0))
+movieplot(ra1, toy.o1, show = 'atypes', adata.show = TRUE,
+ link.col.show = FALSE, link.lty = 0, axes = FALSE,
+ postfn = function(iter) box())
+par(opar)
+
+
+opar <- par(mfrow = c(5, 3), mar = c(0, 0, 0, 0))
+movieplot(ra1, toy.o1, show = 'rwdata', axes = FALSE,
+ postfn = function(iter) box())
+par(opar)
+
+
+
+## ... data set without outliers:
+
+opar <- par(mfrow = c(5, 5), mar = c(0, 0, 0, 0))
+movieplot(a2, toy, show = 'atypes', adata.show = TRUE,
+ link.col.show = FALSE, link.lty = 0, axes = FALSE,
+ postfn = function(iter) box())
+par(opar)
+
+
+par(mfrow = c(5, 5), mar = c(0, 0, 0, 0))
+movieplot(a2, toy, show = 'rwdata', axes = FALSE,
+ postfn = function(iter) box())
+par(opar)
+
+
+
+### Binary bisquare0 reweights function:
+
+set.seed(1235)
+ra2 <- archetypes(toy.o1, 3, family = archetypesFamily('robust',
+ reweightsfn = binary.bisquare0.reweightsfn))
+
+plot(ra2, toy.o1, adata.show = TRUE)
+
+weights(ra2, type = 'reweights')
+
+
+par(mfrow = c(2, 5), mar = c(0, 0, 0, 0))
+movieplot(ra2, toy.o1, show = 'rwdata', axes = FALSE,
+ postfn = function(iter) box())
+par(opar)
+
+
+
+### Compared against leastwise tricube reweightsfn on a
+### data set without outliers:
+
+set.seed(1234)
+a3 <- archetypes(toy, 3, family = archetypesFamily('robust',
+ reweightsfn = leastwise.tricube.reweightsfn))
+plot(a3, toy)
+
+weights(a3, type = 'reweights')
+
+
+par(mfrow = c(5, 4), mar = c(0, 0, 0, 0))
+movieplot(a3, toy, show = 'atypes', adata.show = TRUE,
+ link.col.show = FALSE, link.lty = 0, axes = FALSE,
+ postfn = function(iter) box())
+par(opar)
+
+
+par(mfrow = c(5, 4), mar = c(0, 0, 0, 0))
+movieplot(a3, toy, show = 'rwdata', axes = FALSE,
+ postfn = function(iter) box())
+par(opar)
+
+
+
+### Archetypes' alpine world:
+
+## Panorma view in the case of the original algorithm:
+set.seed(1234)
+a4 <- archetypes(toy.o1, 3)
+
+plot(a4, toy.o1)
+
+panorama(a4, toy.o1)
+panorama(a4, toy, ref.order = 1)
+
+
+## Panorma view in the case of the robust algorithm:
+panorama(a2, toy)
+panorama(a2, toy, ref.order = 1)
+
+par(mfrow = c(3, 1))
+plot(alphas(a4)[, 1])
+plot(alphas(a4)[, 2])
+plot(alphas(a4)[, 3])
+
+
+## A strange idea:
+archetypes.distance.diagplot(a2, toy)
+
+
+
Property changes on: branches/pkg-robust/sandbox/demo2.R
___________________________________________________________________
Name: svn:keywords
+ Date Revision Author URL Id
Name: svn:eol-style
+ native
More information about the Archetypes-commits
mailing list