[Archetypes-commits] r24 - in branches/pkg-robust: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Feb 19 18:26:00 CET 2010
Author: manuel
Date: 2010-02-19 18:26:00 +0100 (Fri, 19 Feb 2010)
New Revision: 24
Added:
branches/pkg-robust/R/archetypes-class-bc.R
branches/pkg-robust/R/archetypes-diagplots.R
Modified:
branches/pkg-robust/DESCRIPTION
branches/pkg-robust/R/archetypes-class.R
branches/pkg-robust/R/archetypes-kit-blocks.R
branches/pkg-robust/R/archetypes-kit.R
branches/pkg-robust/R/archetypes-plot.R
Log:
(1) new interface, (2) weighted + pseudo-robust archetypes
Modified: branches/pkg-robust/DESCRIPTION
===================================================================
--- branches/pkg-robust/DESCRIPTION 2010-02-11 09:26:25 UTC (rev 23)
+++ branches/pkg-robust/DESCRIPTION 2010-02-19 17:26:00 UTC (rev 24)
@@ -1,4 +1,4 @@
-Package: archetypes
+Package: archetypes.dev.robust
Type: Package
Title: Archetypal Analysis
Version: 1.0
Added: branches/pkg-robust/R/archetypes-class-bc.R
===================================================================
--- branches/pkg-robust/R/archetypes-class-bc.R (rev 0)
+++ branches/pkg-robust/R/archetypes-class-bc.R 2010-02-19 17:26:00 UTC (rev 24)
@@ -0,0 +1,214 @@
+.v2 <- function(old, new) {
+ warning(sprintf('Function %s is deprecated; please use %s instead.',
+ sQuote(old), sQuote(new)),
+ call. = FALSE)
+}
+
+#' Archetypes getter.
+#' @param zs An \code{archetypes}-related object.
+#' @param ... Further arguments.
+#' @return Archetypes matrix.
+#' @export
+atypes <- function(zs, ...) {
+ UseMethod('atypes')
+}
+
+#' Archetypes getter.
+#' @param zs An \code{archetypes} object.
+#' @param ... Ignored.
+#' @return Archetypes matrix.
+#' @method atypes archetypes
+#' @S3method atypes archetypes
+atypes.archetypes <- function(zs, ...) {
+ .v2('atypes', 'parameters')
+ return(zs$archetypes)
+}
+
+
+
+#' Number of archetypes getter.
+#' @param zs An \code{archetypes}-related object.
+#' @param ... Further arguments.
+#' @return Number of archetypes.
+#' @export
+ntypes <- function(zs, ...) {
+ UseMethod('ntypes')
+}
+
+#' @S3method ntypes archetypes
+ntypes.archetypes <- function(zs, ...) {
+ return(zs$k)
+}
+
+
+
+#' Residual sum of squares getter.
+#' @param zs An \code{archetypes}-related object.
+#' @param ... Further arguments.
+#' @return Residual sum of squares.
+#' @export
+rss <- function(zs, ...) {
+ UseMethod('rss')
+}
+
+#' Residual sum of squares getter.
+#' @param zs An \code{archetypes} object.
+#' @param ... Ignored.
+#' @return Residual sum of squares.
+#' @method rss archetypes
+#' @S3method rss archetypes
+rss.archetypes <- function(zs, ...) {
+ return(zs$rss)
+}
+
+
+
+#' Archetypes data approximation.
+#' @param zs An \code{archetypes}-related object.
+#' @param ... Further arguments.
+#' @return Approximated data matrix.
+#' @export
+adata <- function(zs, ...) {
+ UseMethod('adata')
+}
+
+#' Archetypes data approximation.
+#' @param zs An \code{archetypes} object.
+#' @param ... Ignored.
+#' @return Approximated data matrix.
+#' @method adata archetypes
+#' @S3method adata archetypes
+adata.archetypes <- function(zs, ...) {
+ .v2('adata', 'fitted')
+ return(t(t(zs$archetypes) %*% t(zs$alphas)))
+}
+
+
+
+#' Alpha getter.
+#' @param zs An \code{archetypes}-related object.
+#' @param ... Further arguments.
+#' @return Alpha matrix.
+#' @export
+alphas <- function(zs, ...) {
+ UseMethod('alphas')
+}
+
+#' Alpha getter.
+#' @param zs An \code{archetypes} object.
+#' @param ... Ignored.
+#' @return Alpha matrix.
+#' @method alphas archetypes
+#' @S3method alphas archetypes
+alphas.archetypes <- function(zs, ...) {
+ .v2('alphas', 'coef')
+ return(zs$alphas)
+}
+
+
+
+#' Beta getter.
+#' @param zs An \code{archetypes}-related object.
+#' @param ... Further arguments.
+#' @return Beta matrix.
+#' @export
+betas <- function(zs, ...) {
+ UseMethod('betas')
+}
+
+#' Beta getter.
+#' @param zs An \code{archetypes} object.
+#' @param ... Ignored.
+#' @return Beta matrix.
+#' @method betas archetypes
+#' @S3method betas archetypes
+betas.archetypes <- function(zs, ...) {
+ .v2('betas', 'coef')
+ return(zs$betas)
+}
+
+
+
+#' Iteration getter.
+#' @param zs An \code{archetypes}-related object.
+#' @param ... Further arguments.
+#' @return Number of iterations.
+#' @export
+iters <- function(zs, ...) {
+ UseMethod('iters')
+}
+
+#' Iteration getter.
+#' @param zs An \code{archetypes} object.
+#' @param ... Ignored.
+#' @return Number of iterations.
+#' @method iters archetypes
+#' @S3method iters archetypes
+iters.archetypes <- function(zs, ...) {
+ return(zs$iters)
+}
+
+
+
+#' Archetypes history getter.
+#' @param zs An \code{archetypes}-related object.
+#' @param ... Further arguments.
+#' @return The \code{archetypes} object of the requested step.
+#' @export
+ahistory <- function(zs, ...) {
+ UseMethod('ahistory')
+}
+
+#' Archetypes history getter.
+#' @param zs An \code{archetypes} object.
+#' @param step The step number to return.
+#' @param ... Ignored.
+#' @return The \code{archetypes} object of the requested step.
+#' @method ahistory archetypes
+#' @S3method ahistory archetypes
+ahistory.archetypes <- function(zs, step, ...) {
+ if ( is.null(zs$history) )
+ stop('No history available')
+
+ if ( step >= 0 )
+ s <- paste('s', step, sep='')
+ else
+ s <- paste('s', nhistory(zs) + step - 1, sep='')
+
+ return(zs$history[[s]][[1]])
+}
+
+
+
+#' Number of history steps getter.
+#' @param zs An \code{archetypes}-related object.
+#' @param ... Further arguments.
+#' @return The number of history steps available.
+#' @export
+nhistory <- function(zs, ...) {
+ UseMethod('nhistory')
+}
+
+#' Archetypes number of history steps getter.
+#' @param zs An \code{archetypes} object.
+#' @param ... Ignored.
+#' @return The number of history steps available.
+#' @method nhistory archetypes
+#' @S3method nhistory archetypes
+nhistory.archetypes <- function(zs, ...) {
+ if ( is.null(zs$history) )
+ stop('No history available')
+
+ return(length(zs$history))
+}
+
+
+#' Kappa getter.
+#' @param z An \code{archetypes} object.
+#' @param ... Ignored.
+#' @return A vector of kappas.
+#' @method kappa archetypes
+#' @S3method kappa archetypes
+kappa.archetypes <- function(z, ...) {
+ return(z$kappas)
+}
Modified: branches/pkg-robust/R/archetypes-class.R
===================================================================
--- branches/pkg-robust/R/archetypes-class.R 2010-02-11 09:26:25 UTC (rev 23)
+++ branches/pkg-robust/R/archetypes-class.R 2010-02-19 17:26:00 UTC (rev 24)
@@ -14,30 +14,41 @@
#' @param kappas The kappas for each system of linear equations.
#' @param betas The data coefficients; a \eqn{p \times n} matrix.
#' @param zas The temporary archetypes.
+#' @param family The archetypes family.
+#' @param residuals The residuals.
+#' @param weights The data weights.
#' @return A list with an element for each parameter and class attribute
#' \code{archetypes}.
#' @seealso \code{\link{archetypes}}, \code{\link{atypes}}, \code{\link{ntypes}},
#' \code{\link{rss}}, \code{\link{adata}}, \code{\link{alphas}},
#' \code{\link{ahistory}}, \code{\link{nhistory}}
#' @export
-as.archetypes <- function(archetypes, k, alphas, rss, iters=NULL, call=NULL,
- history=NULL, kappas=NULL, betas=NULL, zas=NULL) {
-
- return(structure(list(archetypes=archetypes,
- k=k,
- alphas=alphas,
- rss=rss,
- iters=iters,
- kappas=kappas,
- betas=betas,
- zas=zas,
- call=call,
- history=history),
- class='archetypes'))
+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) {
+
+ return(structure(list(archetypes = archetypes,
+ k = k,
+ alphas = alphas,
+ rss = rss,
+ iters = iters,
+ kappas = kappas,
+ betas = betas,
+ zas = zas,
+ call = call,
+ history = history,
+ family = family,
+ residuals = residuals,
+ weights = weights,
+ reweights = reweights),
+ class = c(family$class, 'archetypes')))
}
+setOldClass('archetypes')
+
#' Print method for archetypes object.
#' @param x An \code{archetypes} object.
#' @param full Full information or just convergence and rss information.
@@ -45,219 +56,134 @@
#' @return Undefined.
#' @method print archetypes
#' @S3method print archetypes
-print.archetypes <- function(x, full=TRUE, ...) {
+print.archetypes <- function(x, full = TRUE, ...) {
if ( full ) {
cat('Archetypes object\n\n')
- cat(deparse(x$call), '\n\n')
+ cat(paste(deparse(x$call), collapse = '\n'), '\n\n')
}
-
+
cat('Convergence after', x$iters, 'iterations\n')
- cat('with RSS = ', rss(x), '.\n', sep='')
+ cat('with RSS = ', rss(x), '.\n', sep = '')
}
-#' Archetypes getter.
-#' @param zs An \code{archetypes}-related object.
-#' @param ... Further arguments.
-#' @return Archetypes matrix.
-#' @export
-atypes <- function(zs, ...) {
- UseMethod('atypes')
+#' Return fitted data, i.e. archetypes data approximation.
+#' @param object An \code{archetypes}-related object.
+#' @param ... Ignored.
+#' @return Approximated data matrix.
+#' @method fitted archetypes
+#' @S3method fitted archetypes
+fitted.archetypes <- function(object, ...) {
+ t(t(object$archetypes) %*% t(object$alphas))
}
-#' Archetypes getter.
-#' @param zs An \code{archetypes} object.
+
+
+#' Return fitted archetypes.
+#' @param object An \code{archetypes} object.
#' @param ... Ignored.
#' @return Archetypes matrix.
-#' @method atypes archetypes
-#' @S3method atypes archetypes
-atypes.archetypes <- function(zs, ...) {
- return(zs$archetypes)
+#' @method parameters archetypes
+#' @S3method parameters archetypes
+parameters.archetypes <- function(object, ...) {
+ object$archetypes
}
+#' @importFrom modeltools parameters
+setMethod('parameters', 'archetypes', parameters.archetypes)
-#' Number of archetypes getter.
-#' @param zs An \code{archetypes}-related object.
-#' @param ... Further arguments.
-#' @return Number of archetypes.
-#' @export
-ntypes <- function(zs, ...) {
- UseMethod('ntypes')
+
+#' Return coefficients.
+#' @param object An \code{archetypes} object.
+#' @param type Return alphas or betas.
+#' @param ... Ignored.
+#' @return Coefficient matrix.
+#' @method coef archetypes
+#' @S3method coef archetypes
+coef.archetypes <- function(object, type = c('alphas', 'betas'), ...) {
+ type <- match.arg(type)
+ object[[type]]
}
-#' @S3method ntypes archetypes
-ntypes.archetypes <- function(zs, ...) {
- return(zs$k)
+
+#' Return residuals.
+#' @param object An \code{archetypes} object.
+#' @param ... Ignored.
+#' @return Residuals.
+#' @method residuals archetypes
+#' @S3method residuals archetypes
+residuals.archetypes <- function(object, ...) {
+ object$residuals
}
-#' Residual sum of squares getter.
-#' @param zs An \code{archetypes}-related object.
-#' @param ... Further arguments.
+#' Residual sum of squares.
+#' @param object An \code{archetypes}-related object.
+#' @param ... Ignored.
#' @return Residual sum of squares.
#' @export
-rss <- function(zs, ...) {
+rss <- function(object, ...) {
UseMethod('rss')
}
#' Residual sum of squares getter.
-#' @param zs An \code{archetypes} object.
+#' @param object An \code{archetypes} object.
+#' @param type Return scaled, single or global RSS.
#' @param ... Ignored.
#' @return Residual sum of squares.
#' @method rss archetypes
#' @S3method rss archetypes
-rss.archetypes <- function(zs, ...) {
- return(zs$rss)
-}
+rss.archetypes <- function(object, type = c('scaled', 'single', 'global')) {
+ type <- match.arg(type)
+ resid <- residuals(object)
-
-
-#' Archetypes data approximation.
-#' @param zs An \code{archetypes}-related object.
-#' @param ... Further arguments.
-#' @return Approximated data matrix.
-#' @export
-adata <- function(zs, ...) {
- UseMethod('adata')
+ switch(type,
+ scaled = object$rss,
+ single = apply(resid, 1, object$family$normfn),
+ global = object$family$normfn(resid) / nrow(resid))
}
-#' Archetypes data approximation.
-#' @param zs An \code{archetypes} object.
-#' @param ... Ignored.
-#' @return Approximated data matrix.
-#' @method adata archetypes
-#' @S3method adata archetypes
-adata.archetypes <- function(zs, ...) {
- return(t(t(zs$archetypes) %*% t(zs$alphas)))
-}
-
-#' Alpha getter.
-#' @param zs An \code{archetypes}-related object.
-#' @param ... Further arguments.
-#' @return Alpha matrix.
-#' @export
-alphas <- function(zs, ...) {
- UseMethod('alphas')
+#' Return weights.
+#' @param object An \code{archetypes} object.
+#' @param type Return global weights (weighted archetypes) or
+#' weights calculated during the iterations (robust archetypes).
+#' @return Vector of weights.
+#' @method weights archetypes
+#' @S3method weights archetypes
+weights.archetypes <- function(object, type = c('weights', 'reweights')) {
+ type <- match.arg(type)
+ object[[type]]
}
-#' Alpha getter.
-#' @param zs An \code{archetypes} object.
-#' @param ... Ignored.
-#' @return Alpha matrix.
-#' @method alphas archetypes
-#' @S3method alphas archetypes
-alphas.archetypes <- function(zs, ...) {
- return(zs$alphas)
-}
-
-#' Beta getter.
-#' @param zs An \code{archetypes}-related object.
-#' @param ... Further arguments.
-#' @return Beta matrix.
-#' @export
-betas <- function(zs, ...) {
- UseMethod('betas')
-}
-
-#' Beta getter.
-#' @param zs An \code{archetypes} object.
+#' Predict coefficients or data based on archetypes.
+#' @param object An \code{archetypes} object.
+#' @param type Predict alphas or data.
#' @param ... Ignored.
-#' @return Beta matrix.
-#' @method betas archetypes
-#' @S3method betas archetypes
-betas.archetypes <- function(zs, ...) {
- return(zs$betas)
-}
+#' @return Prediction.
+#' @method predict archetypes
+#' @S3method predict archetypes
+predict.archetypes <- function(object, newdata = NULL,
+ type = c('alphas', 'data'), ...) {
+ type <- match.arg(type)
+ if ( is.null(newdata) )
+ return(switch(type,
+ alphas = coef(object, type = 'alphas'),
+ data = fitted(object)))
+ stop('Not implemented yet.')
-#' Iteration getter.
-#' @param zs An \code{archetypes}-related object.
-#' @param ... Further arguments.
-#' @return Number of iterations.
-#' @export
-iters <- function(zs, ...) {
- UseMethod('iters')
+ ### Something like the following ...
+ #if ( type == 'alphas' )
+ # object$family$alphasfn(NULL, t(object$archetypes), t(newdata))
}
-#' Iteration getter.
-#' @param zs An \code{archetypes} object.
-#' @param ... Ignored.
-#' @return Number of iterations.
-#' @method iters archetypes
-#' @S3method iters archetypes
-iters.archetypes <- function(zs, ...) {
- return(zs$iters)
-}
-
-#' Archetypes history getter.
-#' @param zs An \code{archetypes}-related object.
-#' @param ... Further arguments.
-#' @return The \code{archetypes} object of the requested step.
-#' @export
-ahistory <- function(zs, ...) {
- UseMethod('ahistory')
-}
-
-#' Archetypes history getter.
-#' @param zs An \code{archetypes} object.
-#' @param step The step number to return.
-#' @param ... Ignored.
-#' @return The \code{archetypes} object of the requested step.
-#' @method ahistory archetypes
-#' @S3method ahistory archetypes
-ahistory.archetypes <- function(zs, step, ...) {
- if ( is.null(zs$history) )
- stop('No history available')
-
- if ( step >= 0 )
- s <- paste('s', step, sep='')
- else
- s <- paste('s', nhistory(zs) + step - 1, sep='')
-
- return(zs$history[[s]][[1]])
-}
-
-
-
-#' Number of history steps getter.
-#' @param zs An \code{archetypes}-related object.
-#' @param ... Further arguments.
-#' @return The number of history steps available.
-#' @export
-nhistory <- function(zs, ...) {
- UseMethod('nhistory')
-}
-
-#' Archetypes number of history steps getter.
-#' @param zs An \code{archetypes} object.
-#' @param ... Ignored.
-#' @return The number of history steps available.
-#' @method nhistory archetypes
-#' @S3method nhistory archetypes
-nhistory.archetypes <- function(zs, ...) {
- if ( is.null(zs$history) )
- stop('No history available')
-
- return(length(zs$history))
-}
-
-
-#' Kappa getter.
-#' @param z An \code{archetypes} object.
-#' @param ... Ignored.
-#' @return A vector of kappas.
-#' @method kappa archetypes
-#' @S3method kappa archetypes
-kappa.archetypes <- function(z, ...) {
- return(z$kappas)
-}
Added: branches/pkg-robust/R/archetypes-diagplots.R
===================================================================
--- branches/pkg-robust/R/archetypes-diagplots.R (rev 0)
+++ branches/pkg-robust/R/archetypes-diagplots.R 2010-02-19 17:26:00 UTC (rev 24)
@@ -0,0 +1,51 @@
+
+
+residuals.diagplot <- function(object, ref.order = 1, ...) {
+ y <- residuals(object)
+ x <- seq(length = nrow(y))
+
+ m <- ncol(y)
+
+ ylim <- range(y)
+ ylab <- colnames(y)
+
+ if ( is.null(ref.order) ) {
+ ref.order <- 0
+ ix <- seq(length = nrow(y))
+ }
+ else {
+ ix <- order(y[, ref.order])
+ }
+
+ layout(matrix(seq(length = m), nrow = m, ncol = 1))
+ for ( i in seq(length = m) ) {
+ plot(x, y[ix, i], ylim = ylim,
+ ylab = sprintf('Residuals %s', ylab[i]),
+ xlab = sprintf('Index%s', ifelse(i == ref.order, ' (reference order)', '')),
+ ...)
+ abline(h = 0, lty = 2, col = gray(0.7), ...)
+ }
+
+ invisible(NULL)
+}
+
+
+
+rss.diagplot <- function(object, sort = FALSE, ...) {
+ y <- rss(object, type = 'single')
+ x <- seq(length = length(y))
+
+ if ( sort )
+ y <- sort(y)
+
+ plot(x, y, xlab = 'Index', ylab = 'RSS', ...)
+ abline(h = 0, lty = 2, col = gray(0.7), ...)
+}
+
+
+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', ...)
+}
Modified: branches/pkg-robust/R/archetypes-kit-blocks.R
===================================================================
--- branches/pkg-robust/R/archetypes-kit-blocks.R 2010-02-11 09:26:25 UTC (rev 23)
+++ branches/pkg-robust/R/archetypes-kit-blocks.R 2010-02-19 17:26:00 UTC (rev 24)
@@ -61,11 +61,11 @@
make.dummyfn <- function(huge=200) {
bp.dummyfn <- function(x) {
- y = rbind(x, rep(huge, ncol(x)))
-
+ y = rbind(x, rep(huge, ncol(x)))
+
attr(y, '.Meta') = attr(x, '.Meta')
attr(y, '.Meta')$dummyrow = nrow(y)
-
+
return(y)
}
@@ -79,7 +79,7 @@
#' @return Archetypes zs.
rm.undummyfn <- function(x, zs) {
dr = attr(x, '.Meta')$dummyrow
-
+
return(zs[-dr,])
}
@@ -119,7 +119,7 @@
#' @return The solved linear system.
ginv.zalphasfn <- function(alphas, x) {
require(MASS)
-
+
return(t(ginv(alphas %*% t(alphas)) %*% alphas %*% t(x)))
}
@@ -131,7 +131,7 @@
#' @return The solved linear system.
opt.zalphasfn <- function(alphas, x) {
z <- rnorm(nrow(x)*nrow(alphas))
-
+
fun <- function(z){
z <- matrix(z, ncol=nrow(alphas))
sum( (x - z %*% alphas)^2)
@@ -154,9 +154,9 @@
#' @return Recalculated alpha.
nnls.alphasfn <- function(coefs, C, d) {
require(nnls)
-
+
n = ncol(d)
-
+
for ( j in 1:n )
coefs[,j] = coef(nnls(C, d[,j]))
@@ -175,8 +175,8 @@
nc = ncol(C)
nr = nrow(C)
-
+
s = svd(C, nv=nc)
yint = t(s$u) %*% d
@@ -228,7 +228,7 @@
}
-
+
### Archetypes initialization functions:
#' Init block: generator for random initializtion.
@@ -237,15 +237,15 @@
make.random.initfn <- function(k) {
bp.initfn <- function(x, p) {
-
+
n = ncol(x)
b = matrix(0, nrow=n, ncol=p)
for ( i in 1:p )
b[sample(n, k, replace=FALSE),i] = 1 / k
-
+
a = matrix(1, nrow=p, ncol=n) / p
-
+
return(list(betas=b, alphas=a))
}
@@ -259,12 +259,12 @@
fix.initfn <- function(x, p) {
n = ncol(x)
-
+
b = matrix(0, nrow = n, ncol = p)
b[indizes,] = diag(p)
a = matrix(1, nrow = p, ncol = n) / p
-
+
return(list(betas = b, alphas = a))
}
@@ -273,6 +273,48 @@
+### Weighting functions:
+
+#' Weighting function: move data closer to global center
+#' @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) {
+ if ( is.null(weights) )
+ return(data)
+
+ dr <- attr(data, '.Meta')$dummyrow
+
+ weights <- as.numeric(1 - weights)
+ center <- rowMeans(data[-dr,])
+
+ data[-dr,] <- data[-dr,] + t(weights * t(center - data[-dr,]))
+
+ data
+}
+
+
+
+### Reweighting functions:
+
+#'
+#'
+bisquare0.reweightsfn <- function(resid) {
+ resid <- apply(resid, 2, function(.) sum(abs(.)))
+ resid0 <- resid < sqrt(.Machine$double.eps)
+
+ s <- resid / 6 * median(resid[!resid0])
+
+ ifelse(s < 1, (1 - s^2)^2, 0)
+}
+
+tricube.reweightsfn <- function(resid) {
+ resid <- apply(resid, 2, function(.) sum(abs(.)))
+ ifelse(resid < 1, (1 - resid^3)^3, 0)
+}
+
+
+
### Archetypes family:
#' Archetypes family constructor.
@@ -281,23 +323,59 @@
#' conceptual parts of the algorithm. Currently, only the 'original' family
#' is supported.
#'
-#' @param which The kind of archetypes family; currently ignored.
+#' @param which The kind of archetypes family.
+#' @param ... Exchange blocks predefined by the kind of family.
#' @return A list containing a function for each of the different parts.
#' @seealso \code{\link{archetypes}}
#' @export
-archetypesFamily <- function(which=c('default', 'ginv')) {
- fam <- list(normfn=norm2.normfn,
- scalefn=std.scalefn,
- rescalefn=std.rescalefn,
- dummyfn=make.dummyfn(200),
- undummyfn=rm.undummyfn,
- initfn=make.random.initfn(1),
- alphasfn=nnls.alphasfn,
- betasfn=nnls.betasfn)
+archetypesFamily <- function(which = c('original', 'weighted', 'robust'), ...) {
- fam$zalphasfn <- switch(which[1],
- 'default' = qrsolve.zalphasfn,
- 'ginv' = ginv.zalphasfn)
+ which <- match.arg(which)
+ blocks <- list(...)
- return(fam)
+ family <- do.call(sprintf('.%s.archetypesFamily', which), list())
+ family$which <- which
+ family$which.exchanged <- NULL
+
+ if ( length(blocks) > 0 ) {
+ family$which <- sprintf('%s*', family$which)
+ family$which.exchanged <- names(blocks)
+
+ for ( n in names(blocks) )
+ family[[n]] <- blocks[[n]]
+ }
+
+
+ family
}
+
+.original.archetypesFamily <- function() {
+ list(normfn = norm2.normfn,
+ scalefn = std.scalefn,
+ rescalefn = std.rescalefn,
+ dummyfn = make.dummyfn(200),
+ undummyfn = rm.undummyfn,
+ initfn = make.random.initfn(1),
+ alphasfn = nnls.alphasfn,
+ betasfn = nnls.betasfn,
+ zalphasfn = qrsolve.zalphasfn,
+ weightfn = function(x, weights) x,
+ reweightsfn = function(x, weights) NULL,
+ 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-02-11 09:26:25 UTC (rev 23)
+++ branches/pkg-robust/R/archetypes-kit.R 2010-02-19 17:26:00 UTC (rev 24)
@@ -26,29 +26,36 @@
#' a <- archetypes(toy, 3)
#' @export
#' @note Please see the vignette for a detailed explanation!
-archetypes <- function(data, k, maxIterations=100,
- minImprovement=sqrt(.Machine$double.eps),
- maxKappa=1000, verbose=TRUE, saveHistory=TRUE,
- family=archetypesFamily('default')) {
-
+archetypes <- function(data, k, maxIterations = 100,
+ minImprovement = sqrt(.Machine$double.eps),
+ maxKappa = 1000, verbose = TRUE, saveHistory = TRUE,
+ family = archetypesFamily('original'), weights = NULL) {
+
### Helpers:
mycall <- match.call()
-
+
history <- NULL
- snapshot <- function(name) {
- history[[paste('s', name, sep='')]] <-
- list(archetypes=as.archetypes(t(family$rescalefn(x,
- family$undummyfn(x, zs))), k, alphas=t(alphas),
- betas=t(betas), zas=t(family$rescalefn(x,
- family$undummyfn(x, zas))), rss=rss,
- kappas=kappas))
+ snapshot <- function(i) {
+ history[[sprintf('s%s', i)]] <-
+ 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))
}
+ printIter <- function(i) {
+ cat(i, ': rss = ', formatC(rss, 8, format = 'f'),
+ ', improvement = ', formatC(imp, 8, format = 'f'),
+ '\n', sep = '')
+ }
+
### Data preparation:
- x <- t(data)
- x <- family$scalefn(x)
- x <- family$dummyfn(x)
+ x1 <- t(data)
+ x1 <- family$scalefn(x1)
+ x1 <- family$dummyfn(x1)
+ x0 <- family$weightfn(x1, weights)
+ x <- x0
n <- ncol(x)
m <- nrow(x)
@@ -56,32 +63,40 @@
### Initialization:
init <- family$initfn(x, k)
-
+
betas <- init$betas
alphas <- init$alphas
zs <- x %*% betas
- rss <- family$normfn(zs %*% alphas - x) / n
+ resid <- zs %*% alphas - x
+ rss <- family$normfn(resid) / n
+
zas <- NULL
+ reweights <- NULL
kappas <- c(alphas=kappa(alphas), betas=kappa(betas),
zas=-Inf, zs=kappa(zs))
isIll <- c(kappas) > maxKappa
errormsg <- NULL
-
+
if ( saveHistory ) {
history <- new.env(parent=emptyenv())
snapshot(0)
}
-
+
### Main loop:
i <- 1
imp <- +Inf
tryCatch(while ( (i <= maxIterations) & (imp >= minImprovement) ) {
-
+
+ ## Reweight data:
+ reweights <- family$reweightsfn(resid)
+ x <- family$weightfn(x0, reweights)
+
+
## Alpha's:
alphas <- family$alphasfn(alphas, zs, x)
zas <- family$zalphasfn(alphas, x)
@@ -89,39 +104,40 @@
kappas[c('alphas', 'zas')] <- c(kappa(alphas), kappa(zas))
-
+
## Beta's:
betas <- family$betasfn(betas, x, zas)
zs <- x %*% betas
kappas[c('betas', 'zs')] <- c(kappa(betas), kappa(zs))
-
- ## RSS and improvement:
- rss2 <- family$normfn(zs %*% alphas - x) / n
-
+
+ ## Residuals, RSS and improvement:
+ alphas0 <- family$alphasfn(alphas, zs, x0)
+
+ resid <- zs %*% alphas0 - x0
+ rss2 <- family$normfn(resid) / n
+
imp <- rss - rss2
rss <- rss2
- kappas <- c(alphas=kappa(alphas), betas=kappa(betas),
- zas=kappa(zas), zs=kappa(zs))
+
+ ## Loop Zeugs:
+ kappas <- c(alphas = kappa(alphas), betas = kappa(betas),
+ zas = kappa(zas), zs = kappa(zs))
isIll <- isIll & (kappas > maxKappa)
-
- ## Loop Zeugs:
if ( verbose )
- cat(i, ': rss = ', formatC(rss, 8, format='f'),
- ', improvement = ', formatC(imp, 8, format='f'),
- '\n', sep = '')
-
+ printIter(i)
+
if ( saveHistory )
snapshot(i)
-
+
i <- i + 1
},
error=function(e) errormsg <<- e)
-
+
### Check illness:
if ( !is.null(errormsg) ) {
warning('k=', k, ': ', errormsg)
@@ -133,14 +149,27 @@
warning('k=', k, ': ', paste(names(isIll)[isIll], collapse=', '),
' > maxKappa', sep='')
-
- ### Rescale archetypes:
+
+ ### Rescale archetypes, etc.:
+ if ( !is.null(weights) || !is.null(reweights) ) {
+ alphas <- family$alphasfn(alphas, zs, x1)
+ betas <- family$betasfn(betas, x1, zs)
+ }
+
zs <- family$undummyfn(x, zs)
zs <- family$rescalefn(x, zs)
- zs <- t(zs)
-
- return(as.archetypes(zs, k, t(alphas), rss, iters=(i-1),
- call=mycall, history=history, kappas=kappas,
- betas=t(betas)))
+
+ ### Recalculate residuals, etc. for original data:
+ resid <- zs %*% alphas - t(data)
+
+
+ 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))
}
+
+
+
Modified: branches/pkg-robust/R/archetypes-plot.R
===================================================================
--- branches/pkg-robust/R/archetypes-plot.R 2010-02-11 09:26:25 UTC (rev 23)
+++ branches/pkg-robust/R/archetypes-plot.R 2010-02-19 17:26:00 UTC (rev 24)
@@ -7,13 +7,13 @@
a <- rbind(atypes(zs), atypes(zs)[1,])
xc <- a[,1]; xm <- mean(xc)
yc <- a[,2]; ym <- mean(yc)
-
+
real <- xc - xm
imag <- yc - ym
angle <- atan2(imag, real)
-
+
index <- order(angle)
-
+
return(a[c(index, index[1]),])
}
@@ -46,39 +46,73 @@
#' @export
#' @noRd
plot.archetypes <- function(x, y,
- data.col=gray(0.7), data.pch=19,
- atypes.col=2, atypes.pch=19,
- ahull.show=TRUE, ahull.col=atypes.col,
- chull=NULL, chull.col=1, chull.pch=19,
- adata.show=FALSE, adata.col=3, adata.pch=13,
- link.col=data.col, ...) {
+ data.col = gray(0.7), data.pch = 19, data.bg = NULL,
+ atypes.col = 2, atypes.pch = 19,
+ ahull.show = TRUE, ahull.col = atypes.col,
+ chull = NULL, chull.col = 1, chull.pch = 19,
+ adata.show = FALSE, adata.col = 3, adata.pch = 13,
+ link.col = data.col, link.lty = 1, ...) {
zs <- x; data <- y;
- plot(data, col=data.col, pch=data.pch, ...)
- points(atypes(zs), col=atypes.col, pch=atypes.pch, ...)
+ plot(data, col = data.col, pch = data.pch, bg = data.bg, ...)
+ points(atypes(zs), col = atypes.col, pch = atypes.pch, ...)
if ( !is.null(chull) ) {
- points(data[chull,], col=chull.col, pch=chull.pch, ...)
- lines(data[c(chull, chull[1]),], col=chull.col, ...)
+ points(data[chull,], col = chull.col, pch = chull.pch, ...)
+ lines(data[c(chull, chull[1]),], col = chull.col, ...)
}
if ( ahull.show )
- lines(ahull(zs), col=ahull.col)
+ lines(ahull(zs), col = ahull.col)
if ( adata.show ) {
### Based on an idea of Bernard Pailthorpe.
adata <- adata(zs)
-
- points(adata, col=adata.col, pch=adata.pch, ...)
+ link.col <- rep(link.col, length = nrow(adata))
+ link.lty <- rep(link.lty, length = nrow(adata))
+
+ points(adata, col = adata.col, pch = adata.pch, ...)
for ( i in seq_len(nrow(data)) )
- lines(rbind(data[i,], adata[i,]), col=link.col, ...)
+ lines(rbind(data[i,], adata[i,]), col = link.col[i],
+ lty = link.lty[i], ...)
}
}
+#' Plot weighted archetypes.
+plot.weightedArchetypes <- function(x, y,
+ adata.show = FALSE, data.col = 1, data.pch = 21,
+ data.bg = gray, link.col.show = TRUE,
+ weights.type = 'weights', ...) {
+ if ( adata.show ) {
+ w <- 1 - weights(x, type = weights.type)
+
+ if ( link.col.show ) {
+ link.col <- ifelse(w == 1, 1, data.bg(w))
+ link.lty <- ifelse(w == 1, 2, 1)
+ }
+
+ 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)
+ }
+ else {
+ plot.archetypes(x, y, ...)
+ }
+}
+
+
+
+#' Plot robust archetypes.
+plot.robustArchetypes <- function(x, y, ...) {
+ plot.weightedArchetypes(x, y, weights.type = 'reweights', ...)
+}
+
+
+
#' Plot of data and stepArchetypes.
#' @param x An \code{\link{stepArchetypes}} object.
#' @param y A matrix or data frame.
@@ -97,11 +131,11 @@
data.col=gray(0.7), data.pch=19,
atypes.col=(seq_len(length(x) * length(x[[1]]))+1),
atypes.pch=19, ahull.show=TRUE, ahull.col=atypes.col, ...) {
-
+
zs <- x; data <- y;
-
+
flatzs <- unlist(zs, recursive=FALSE)
-
+
plot(data, col=data.col, pch=data.pch, ...)
for ( i in seq_along(flatzs) ) {
a <- flatzs[[i]]
More information about the Archetypes-commits
mailing list