Clústering

No hay nada tan corrosivo para la fe en el clústering que probar una y otra vez k-medias (por ejemplo) sobre los mismos datos y ver cómo los resultados cambian drásticamente de ejecución en ejecución.

Pero eso viene a ser, esencialmente, lo que hay detrás del consensus clústering (CC), una técnica que puede ser usada, entre otros fines, para determinar el número óptimo de grupos.

La idea fundamental de la cosa es que observaciones que merezcan ser agrupadas juntas lo serán muy frecuentemente aunque cambien ligeramente las condiciones iniciales (por ejemplo, se tome una submuestra de los datos o cambien las condiciones iniciales de k-medias, por ejemplo). Si uno altera esas condiciones iniciales repetidas veces puede contar la proporción de las veces que las observaciones i y j fueron emparejadas juntas y crear la correspondiente matriz (simétrica, para más señas) $latex C(i,j)$.

De vez en cuando pruebo paquetes promisorios. No es infrecuente que cosas que he intentado hace años, algún ejemplo más o menos sencillo que he publicado aquí, acabe convirtiéndose en la piedra angular de algo facturable. Incluso de algo facturable por mí.

geozoning podía haber sido uno de esos. La promesa del paquete es que puede ayudarte a segmentar regiones del espacio de acuerdo con alguna variable, una especie de clústering para información de tipo espacial.

El otro día pasé por uno de esos sitios en los que exponen en las paredes obras de artistas medianos con el precio debajo. Me quedé mirando una muy… concéntrica porque me recordaba a lo que nos regala a menudo Antonio Chinchón. Pregunté de qué trataba la cosa y tuvieron la paciencia de explicármelo: al lado había una foto enorme y, se conoce, las cosas concéntricas eran una reordenación de los píxels de la primera. Una especie de tortilla de patata deconstruida a lo Adriá, pero con fotos.

En mi clase de aprendizaje no supervisado en el máster de ciencia de datos de la U-TAD mostré un vídeo en el que se ilustraba el funcionamiento del algoritmo de las k-medias.

Una alumna encontró un recurso mucho mejor.

Que trae, además, como bonus, una ilustración del funcionamiento de DBSCAN (véase también esto).

Los viejos del lugar recordarán esto, donde agrupo trayectorias usando k-medias a pelo.

El paquete kmlShape usa la distancia de Fréchet para hacer algo parecido: buscar trayectorias geométricamente similares.

El código es

    library(kmlShape)
    library(tseries)
    library(zoo)
    library(XML)
    library(reshape)
    library(ggplot2)

    foo  <- function(
      simbolo, final = Sys.time(),
      profundidad = 30 * 24 * 3600) {
      precios <- get.hist.quote(
        instrument= simbolo,
        start = final - profundidad,
        end = final, quote=c("AdjClose"),
        provider="yahoo", origin="1970-01-01",
        compression="d", retclass="zoo")
      colnames(precios) <- simbolo
      return(precios)
    }

    # lista de símbolos del ibex

    tmp <- readHTMLTable("http://finance.yahoo.com/q/cp?s=%5EIBEX+Components")[[5]]
    tmp <- as.character(tmp$V1[-(1:6)])

    ibex <- do.call(merge,
      sapply(simbolos, foo, simplify = F))

    ibex.scaled <- data.frame(t(scale(ibex)))
    tmp <- cldsWide(ibex.scaled)

    res <- kmlShape(tmp, 4, toPlot = "none")

    tmp <- data.frame(
      id = rownames(ibex.scaled),
      cluster = res@clusters, ibex.scaled)

    tmp <- melt(tmp, id.vars = c("id", "cluster"))
    tmp$fecha <- as.Date(tmp$variable, "X%Y.%m.%d")

    ggplot(tmp, aes(x=fecha, y=value, group=id)) +
      geom_line() + facet_wrap(~cluster)

y el resultado,

El clústering (o análisis de conglomerados, o como se le quiera llamar) es un atajo intelectual. En lugar de estudiar individuos (no necesariamente personas), estos se agrupan de manera más o menos cuestionable, se eligen representantes en cada uno de ellos, cuyas características se atribuyen a continuación a todos sus miembros.

No puedo evitar escribir párrafos como el anterior sin que me venga a la nariz ese olor a naftalina de cuando era crío y abría los armarios de mi abuela.

Ha sido en latitudes otras que las habituales que he aprendido y leído (mas no probado) sobre DBSCAN. Se conoce que es un nuevo (aunque ya tiene sus añitos: algo así como 20) método de clústering.

Por un lado, se agradecen las novedades.

Por el otro, tengo cierta aversión a las cosas que proceden de los congresos de Knowledge Discovery and Data Mining, que es donde fue publicado el algoritmo.

En esencia, funciona así: se fijan dos parámetros, e y n. Un punto es central si a distancia e o menor tiene, al menos, otros n puntos. Los clústers los conforman:

Especulando sobre la diferencia en la práctica entre distintas métricas ($latex l_1$, $latex l_2$, $latex l_\infty$, etc.), construi una serie de diagramas de Voronoi usado métricas arbitrarias.

En la Wikipedia se comparan gráficamente $latex l_1$, $latex l_2$ (o euclídea y Manhattan). Mi código,

library(data.table)
library(reshape2)
library(grid)

n <- 20
dim.image <- 1000
puntos <- data.frame(id = 1:n,
                      x0 = runif(n) * dim.image,
                      y0 = runif(n) * dim.image)
colores <- rainbow(n)

voronoi <- function(p){
  tmp <- data.table(expand.grid(
      x = 1:dim.image,
      y = 1:dim.image, id = 1:n), key = "id")
  tmp <- merge(tmp, puntos, by = "id")

  distancia <- function(a, b, c, d, p)
    (abs(a-c)^p + abs(b-d)^p)^(1/p)

  tmp$distancia <- distancia(tmp$x,
    tmp$y, tmp$x0, tmp$y0, p)
  tmp[, rank := rank(distancia, ties = "random"),
    by = c("x", "y")]

  rejilla <- tmp[tmp$rank == 1,]
  rejilla$x0 <- rejilla$y0 <-
    rejilla$distancia <- rejilla$rank <- NULL

  rejilla$color <- colores[rejilla$id]

  imagen <- as.matrix(dcast(rejilla, x ~ y, value.var = "color")[,-1])

  grid.raster(imagen)
}

permite usar más en función del parámetro p.

A cuento de la sesión práctica de modelos no supervisados que impartiré este sábado y que estoy preparando justo ahora, traigo a la atención de mis lectores una imagen que el asunto me sugiere:

La fuente también vale la pena. Aunque habla de otra cosa.

Nota: releyendo la entrada antes de publicarla definitivamente, me doy cuenta de que igual estoy siendo excesivamente sutil.

Te escribo desde Zúrich. Tenía, o eso creía yo, un vuelo a Madrid esta tarde. Pero al llegar al aeropuerto, no constaba en la lista de viajeros. El motivo figura en una subsubsubpágina de vuestro portal: si un viajero pierde el vuelo de ida, como me ocurrió el lunes pasado, se cancela automáticamente el de vuelta.

Como tú también has viajado mucho, podrás imaginar lo que ha sucedido en el aeropuerto: paseos de mostrador en mostrador hasta recibir una oferta rayana en impuesto revolucionario: me han ofrecido uno de esos asientos vacíos por un precio que, lo siento, no estoy dispuesto a pagar. He vuelto a casa y he comprado un billete a la competencia: volaré desde Basilea, la ciudad de nuestros bienamados Bernoullis, el sábado a las 6:10 AM. Viajaré en último tren del viernes y mataré el rato en el aeropuerto como cuando estudiante.

Guardo desde hace un tiempo el enlace al paquete blockcluster de R que igual puede ser del interés de alguno de mis lectores.

No lo he probado pero sospecho que cualquier día me puede sacar de un apuro. Implementa lo que dice, el coclústering, concepto que se explica mejor, como el efecto de las dietas milagrosas, con la foto del antes y el después:

Esto es: la entrada es una matriz y la salida es una matriz reorganizada tanto en sus filas como en sus columnas en la que se han detectado bloques homogéneos.

Estar suscrito a las actualizaciones de CRAN le permite a uno estar al tanto de las novedades de R de otra manera. De vez en cuando uno encuentra pequeños paquetes que le solucionan un problema puntual. Mucho más frecuentemente, la verdad, uno se topa con aplicaciones muy específicas en áreas que le resultan remotas.

Pero uno no espera nunca tropiezar con paquetes que no sabe si clasificar como una broma, una ironía bromas o como algo mucho peor: la constatación de una triste realidad. Es el caso de ykmeans.

Inspirado en esto aunque con la intención de mejorar el horrible código adjunto, escribí el otro día esto:

library("biOps")
library("cluster")

# leo una foto usando readJpeg de biOps
# el objeto devuelto es un array mxnx3 dimensional
# la última dimensión es el rgb de cada pixel

tmp <- tempfile()
download.file("http://blog.guiasenior.com/images/Retrato_Garber.jpg", tmp)
x <- readJpeg(tmp)

# si quieres mostrar la foto como un gráfico...
#plot(x)

# convertimos el array 3D nxmx3 en uno 2D (nm)x3
# luego buscamos 5 clústers
# esencialmente, buscamos 7 "píxels representativos"
d <- dim(x)
clarax <- clara(array(x, dim = c(d[1] * d[2], d[3])), 7)

# reemplazamos cada rgb de cada cluster por su
# "píxel representativo" (medioide) correspondiente
rgb.clusters <- clarax$medoids[clarax$cluster,]

# convertimos la matriz resultante en un array 3D
# (invirtiendo la transformación anterior)
# y representamos gráficamente
plot(imagedata(array(rgb.clusters, dim = d)))

Obviamente, podéis cambiar la foto y hacer variar el número de clústers. Pero conviene recordar que: