[Archetypes-commits] r20 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Feb 8 14:05:05 CET 2010
Author: manuel
Date: 2010-02-08 14:05:05 +0100 (Mon, 08 Feb 2010)
New Revision: 20
Modified:
pkg/R/archetypes-class.R
pkg/R/archetypes-kit-blocks.R
pkg/R/archetypes-kit.R
pkg/R/archetypes-pcplot.R
pkg/R/archetypes-plot.R
Log:
weighted archetypes.
Modified: pkg/R/archetypes-class.R
===================================================================
--- pkg/R/archetypes-class.R 2009-04-27 08:27:34 UTC (rev 19)
+++ pkg/R/archetypes-class.R 2010-02-08 13:05:05 UTC (rev 20)
@@ -14,6 +14,7 @@
#' @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 weights Vector of data weights within \eqn{[0, 1]}.
#' @return A list with an element for each parameter and class attribute
#' \code{archetypes}.
#' @seealso \code{\link{archetypes}}, \code{\link{atypes}}, \code{\link{ntypes}},
@@ -21,19 +22,26 @@
#' \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'))
+ history=NULL, kappas=NULL, betas=NULL, zas=NULL,
+ weights=NULL) {
+
+ a <- structure(list(archetypes=archetypes,
+ k=k,
+ alphas=alphas,
+ rss=rss,
+ iters=iters,
+ kappas=kappas,
+ betas=betas,
+ zas=zas,
+ call=call,
+ history=history,
+ weights=weights),
+ class='archetypes')
+
+ if ( !is.null(weights) )
+ class(a) <- c('weightedArchetypes', class(a))
+
+ a
}
@@ -50,7 +58,7 @@
cat('Archetypes object\n\n')
cat(deparse(x$call), '\n\n')
}
-
+
cat('Convergence after', x$iters, 'iterations\n')
cat('with RSS = ', rss(x), '.\n', sep='')
}
@@ -223,7 +231,7 @@
s <- paste('s', step, sep='')
else
s <- paste('s', nhistory(zs) + step - 1, sep='')
-
+
return(zs$history[[s]][[1]])
}
Modified: pkg/R/archetypes-kit-blocks.R
===================================================================
--- pkg/R/archetypes-kit-blocks.R 2009-04-27 08:27:34 UTC (rev 19)
+++ pkg/R/archetypes-kit-blocks.R 2010-02-08 13:05:05 UTC (rev 20)
@@ -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,23 @@
+### 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)
+
+ weights <- as.numeric(1 - weights)
+ center <- rowMeans(data)
+ data + t(weights * t(center - data))
+}
+
+
+
### Archetypes family:
#' Archetypes family constructor.
@@ -293,7 +310,8 @@
undummyfn=rm.undummyfn,
initfn=make.random.initfn(1),
alphasfn=nnls.alphasfn,
- betasfn=nnls.betasfn)
+ betasfn=nnls.betasfn,
+ weightfn=center.weightfn)
fam$zalphasfn <- switch(which[1],
'default' = qrsolve.zalphasfn,
Modified: pkg/R/archetypes-kit.R
===================================================================
--- pkg/R/archetypes-kit.R 2009-04-27 08:27:34 UTC (rev 19)
+++ pkg/R/archetypes-kit.R 2010-02-08 13:05:05 UTC (rev 20)
@@ -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 Vector of data weights within \eqn{[0, 1]}.
#' @param maxIterations The maximum number of iterations.
#' @param minImprovement The minimal value of improvement between two
#' iterations.
@@ -26,14 +27,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('default')) {
-
+
### Helpers:
mycall <- match.call()
-
+
history <- NULL
snapshot <- function(name) {
history[[paste('s', name, sep='')]] <-
@@ -41,14 +42,15 @@
family$undummyfn(x, zs))), k, alphas=t(alphas),
betas=t(betas), zas=t(family$rescalefn(x,
family$undummyfn(x, zas))), rss=rss,
- kappas=kappas))
+ kappas=kappas, weights=weights))
}
### Data preparation:
- x <- t(data)
- x <- family$scalefn(x)
- x <- family$dummyfn(x)
+ data <- t(data)
+ data <- family$scalefn(data)
+ data <- family$dummyfn(data)
+ x <- family$weightfn(data, weights)
n <- ncol(x)
m <- nrow(x)
@@ -56,7 +58,7 @@
### Initialization:
init <- family$initfn(x, k)
-
+
betas <- init$betas
alphas <- init$alphas
@@ -69,19 +71,19 @@
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) ) {
-
+
## Alpha's:
alphas <- family$alphasfn(alphas, zs, x)
zas <- family$zalphasfn(alphas, x)
@@ -89,39 +91,39 @@
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
-
+
imp <- rss - rss2
rss <- rss2
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 = '')
-
+
if ( saveHistory )
snapshot(i)
-
+
i <- i + 1
},
error=function(e) errormsg <<- e)
-
+
### Check illness:
if ( !is.null(errormsg) ) {
warning('k=', k, ': ', errormsg)
@@ -133,14 +135,20 @@
warning('k=', k, ': ', paste(names(isIll)[isIll], collapse=', '),
' > maxKappa', sep='')
-
+
### Rescale archetypes:
+ if ( !is.null(weights) ) {
+ alphas <- family$alphasfn(alphas, zs, data)
+ betas <- family$betasfn(betas, data, 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)))
+ betas=t(betas), weights=weights))
}
+
Modified: pkg/R/archetypes-pcplot.R
===================================================================
--- pkg/R/archetypes-pcplot.R 2009-04-27 08:27:34 UTC (rev 19)
+++ pkg/R/archetypes-pcplot.R 2010-02-08 13:05:05 UTC (rev 20)
@@ -7,6 +7,7 @@
#' @param data A matrix or data frame.
#' @param data.col Color of data lines.
#' @param data.lwd Width of data lines.
+#' @param data.lty Type of data lines.
#' @param atypes.col Color of archetypes lines.
#' @param atypes.lwd Width of archetypes lines.
#' @param atypes.lty Type of archetypes lines.
@@ -19,17 +20,36 @@
#' @return Undefined.
#' @method pcplot archetypes
#' @S3method pcplot archetypes
-pcplot.archetypes <- function(x, data, data.col=gray(0.7), data.lwd=1,
+pcplot.archetypes <- function(x, data, data.col=gray(0.7), data.lwd=1, data.lty=1,
atypes.col=2, atypes.lwd=2, atypes.lty=1,
chull=NULL, chull.col=1, chull.lwd=2, chull.lty=1, ...) {
- pcplot(data, col=data.col, lwd=data.lwd, ...)
+ pcplot(data, col=data.col, lwd=data.lwd, lty=data.lty, ...)
if ( !is.null(chull) )
lines.pcplot(data[chull,], data,
col=chull.col, lwd=chull.lwd, lty=chull.lty, ...)
-
+
lines.pcplot(atypes(x), data,
col=atypes.col, lwd=atypes.lwd, lty=atypes.lty, ...)
}
+
+
+#' Parallel coordinates of weighted data and archetypes.
+#' @param x An \code{\link{archetypes}} object.
+#' @param data A matrix or data frame.
+#' @param data.col Function to calculate weighted data lines color.
+#' @param ... Passed to \code{\link{pcplot.archetypes}}.
+#' @return Undefined.
+#' @method pcplot weightedArchetypes
+#' @S3method pcplot weightedArchetypes
+#' @noRd
+pcplot.weightedArchetypes <- function(x, data,
+ data.col=function(x) gray(1 - x), ...) {
+
+ col <- data.col(x$weights)
+ lty <- ifelse(x$weights == 0, 2, 1)
+
+ pcplot.archetypes(x, data, data.col=col, data.lty=lty, ...)
+}
Modified: pkg/R/archetypes-plot.R
===================================================================
--- pkg/R/archetypes-plot.R 2009-04-27 08:27:34 UTC (rev 19)
+++ pkg/R/archetypes-plot.R 2010-02-08 13:05:05 UTC (rev 20)
@@ -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]),])
}
@@ -24,6 +24,7 @@
#' @param y A matrix or data frame.
#' @param data.col Color of data points.
#' @param data.pch Type of data points.
+#' @param data.bg Fill color for data points.
#' @param atypes.col Color of archetypes points.
#' @param atypes.pch Type of archetypes points.
#' @param ahull.show Show approximated convex hull.
@@ -46,7 +47,7 @@
#' @export
#' @noRd
plot.archetypes <- function(x, y,
- data.col=gray(0.7), data.pch=19,
+ 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,
@@ -55,7 +56,7 @@
zs <- x; data <- y;
- plot(data, col=data.col, pch=data.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) ) {
@@ -70,15 +71,41 @@
if ( adata.show ) {
### Based on an idea of Bernard Pailthorpe.
adata <- adata(zs)
-
+ link.col <- rep(link.col, length.out=nrow(adata))
+
+ for ( i in seq_len(nrow(data)) )
+ lines(rbind(data[i,], adata[i,]), col=link.col[i], ...)
+
points(adata, col=adata.col, pch=adata.pch, ...)
- for ( i in seq_len(nrow(data)) )
- lines(rbind(data[i,], adata[i,]), col=link.col, ...)
}
+
+
+ invisible(NULL)
}
+#' Plot of weighted data and archetypes.
+#' @param x An \code{\link{archetypes}} object.
+#' @param y A matrix or data frame.
+#' @param data.col Color of data points.
+#' @param data.pch Type of data points.
+#' @param data.bg Function to calculate weighted data point color.
+#' @param ... Passed to the underlying \code{\link{plot.archetypes}} functions.
+#' @return Undefined.
+#' @method plot weightedArchetypes
+#' @export
+#' @noRd
+plot.weightedArchetypes <- function(x, y,
+ data.col=1, data.pch=21,
+ data.bg=function(x) gray(1 - x), ...) {
+
+ plot.archetypes(x, y, data.pch = data.pch, data.col=data.col,
+ data.bg=data.bg(x$weights), ...)
+}
+
+
+
#' Plot of data and stepArchetypes.
#' @param x An \code{\link{stepArchetypes}} object.
#' @param y A matrix or data frame.
@@ -97,11 +124,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]]
@@ -113,3 +140,4 @@
}
+
More information about the Archetypes-commits
mailing list