Ciencia De Datos

En algunas de las últimas charlas (de ML) a las que he asistido se han enumerado recetas con las que tratar de resolver distintos problemas. Pero no han explicado cuándo ni por qué es conveniente aplicarlas. Incluso cuando se han presentado dos y hasta tres recetas para el mismo problema.

Me consta que parte de la audiencia quedó desconcertada y falta de algo más. ¿Tal vez una receta para aplicar recetas? ¿De una metarreceta?

De los últimos foros de científicos de datos a los que he asistido, de las últimas conversaciones con científicos de datos que he mantenido, he salido con una gran duda: ¿soy yo el que tiende a juntarse con ellos o es que hay una plaga de talibanes del RMSE es Dios y XGBoost, su profeta?

Lejos está ese lema simplificador de los principios que me mueven a escribir estas páginas. Por lo que, anuncio, estoy arrejuntando razones y papelotes con los que tratar de arrancar un movimiento herético.

Ese es el nombre agramatical de una nueva consejería pergeñada por una red neuronal recurrente que he ajustado usando un año de BOEs.

El código, adaptado de aquí y sustancialmente mejorado, es

library(mxnet)

batch.size     <- 32
seq.len        <- 64
num.hidden     <- 128
num.embed      <- 8
num.lstm.layer <- 1
num.round      <- 1
learning.rate  <- 0.1
wd             <- 0.00001
clip_gradient  <- 1
update.period  <- 1

make.data <- function(dir.boe, seq.len = 32,
  max.vocab=10000, dic = NULL) {
  text <- lapply(dir(dir.boe), readLines)
  text <- lapply(text, paste, collapse = "\n")
  text <- paste(text, collapse = "\n")

  char.lst <- strsplit(text, '')[[1]]
  chars <- unique(char.lst)

  num.seq  <- floor(length(char.lst) / seq.len)
  char.lst <- char.lst[1:(num.seq * seq.len)]
  data <- matrix(match(char.lst, chars) - 1, seq.len, num.seq)

  dic <- as.list(1:length(chars))
  names(dic) <- chars

  lookup.table <- as.list(chars)

  return (list(data = data, dic = dic,
    lookup.table = lookup.table))
}


ret <- make.data(".", seq.len=seq.len)

X   <- ret$data
dic <- ret$dic
lookup.table <- ret$lookup.table

vocab <- length(dic)

train.val.fraction <- 0.9
train.cols <- floor(ncol(X) * train.val.fraction)

drop.tail <- function(x, batch.size) {
  nstep <- floor(ncol(x) / batch.size)
  x[, 1:(nstep * batch.size)]
}

get.label <- function(X)
  matrix(c(X[-1], X[1]), nrow(X), ncol(X))

X.train.data   <- X[, 1:train.cols]
X.train.data   <- drop.tail(X.train.data, batch.size)
X.train.label  <- get.label(X.train.data)
X.train        <- list(data=X.train.data, label=X.train.label)

X.val.data     <- X[, -(1:train.cols)]
X.val.data     <- drop.tail(X.val.data, batch.size)
X.val.label    <- get.label(X.val.data)
X.val          <- list(data=X.val.data, label=X.val.label)


model <- mx.lstm(X.train, X.val,
    ctx=mx.cpu(),
    num.round=num.round,
    update.period=update.period,
    num.lstm.layer=num.lstm.layer,
    seq.len=seq.len,
    num.hidden=num.hidden,
    num.embed=num.embed,
    num.label=vocab,
    batch.size=batch.size,
    input.size=vocab,
    initializer=mx.init.uniform(0.1),
    learning.rate=learning.rate,
    wd=wd,
    clip_gradient=clip_gradient)


get.sample <- function(n, start = "<", random.sample = TRUE){

  make.output <- function(prob, sample = FALSE) {
    prob <- as.numeric(as.array(prob))
    if (!sample)
      return(which.max(as.array(prob)))
    sample(1:length(prob), 1, prob = prob^2)
  }

  infer.model <- mx.lstm.inference(
      num.lstm.layer=num.lstm.layer,
      input.size=vocab,
      num.hidden=num.hidden,
      num.embed=num.embed,
      num.label=vocab,
      arg.params=model$arg.params,
      ctx=mx.cpu())

  out <- start
  last.id <- dic[[start]]

  for (i in 1:(n-1)) {
    ret <- mx.lstm.forward(infer.model, last.id - 1, FALSE)
    infer.model <- ret$model
    last.id <- make.output(ret$prob, random.sample)
    out <- paste0(out, lookup.table[[last.id]])
  }
  out
}

cat(get.sample(1000, start = "A", random.sample = T))

Lo anterior genera cosas tales como:

El otro día asistí a la enésima confusión sobre k-medias y k-vecinos. Que lo es, más en general, sobre el clústering contra modelos locales de la clase que sean, desde k-vecinos hasta el filtrado colaborativo. Veamos si esta comparación que traigo hoy a mis páginas contribuye a erradicar dicha confusión.

k-medias es como las elecciones. Hace poco tuvimos unas en España. Alguien decidió (aproximadamente) que k = 4 y nos pidió, a nosotros, punticos del espacio, identificar el centroide más próximo a nosotros para que lo votásemos. Pues eso, la misma frustración que muchos dizque sintieron teniendo que elegir entre partidos/centroides subjetivamente igual de alejados de los intereses de uno es la que sienten nuestros punticos cuando los procrusteamos para asociarlos al totum revolutum de los clientes estrella, etc.

En esto.

Durante la charla de Carlos Ortega del pasado jueves sobre el paquete caret y sus concomitancias, se planteó el asunto de la optimización de los parámetros de un modelo usando rejillas (grids) de búsqueda.

Cuando un determinado algoritmo depende de, p.e., cuatro parámetros, se puede definir una rejilla como en

gbmGrid <-  expand.grid(interaction.depth = c(1, 5, 9),
      n.trees = (1:30)*50,
      shrinkage = 0.1,
      n.minobsinnode = 20)

y caret se encarga de ajustar el modelo bajo todas esas combinaciones de parámetros (90 en el ejemplo) para ver cuál de ellas es, con las debidas salvedades, óptima.

Imaginemos un banco que construye modelos para determinar si se concede o no un crédito. Este banco tiene varias opciones para crear el modelo. Sin embargo, en algunos países el regulador exige que el banco pueda explicar el motivo de la denegación de un crédito cuando un cliente lo solicite.

Esa restricción impediría potencialmente usar modelos de caja negra como el que construyo a continuación:

library(randomForest)

raw <- read.table("http://archive.ics.uci.edu/ml/machine-learning-databases/credit-screening/crx.data",
    sep = ",", na.strings = "?")

dat <- raw
dat$V14 <- dat$V6 <- NULL    # me da igual
dat <- na.omit(dat)          # ídem

modelo <- randomForest(V16 ~ ., data = dat)

Fijémonos en el sujeto 100, a quien se le deniega el crédito (suponiendo, ¡mal hecho!, que el punto de corte de la probabilidad para concederlo es el 50%), y la variable $V8$. Podemos ver cuál sería el score del cliente modificando esa variable entre su valor real y el máximo del rango dejando las demás tal cual:

Es

$$ \sum_i \Phi(y_i, f_1(x_i)) > \sum_i \Phi(y_i, f_1(x_i) - \lambda \nabla \Phi(y_i, f_1(x_i)) \sim$$ $$ \sim \sum_i \Phi(y_i, f_1(x_i) - \lambda f_2(x_i))$$

Por supuesto, el lector se preguntará muchas cosas, entre las que destaco:

• ¿Qué representa cada uno de los elementos que aparecen en la línea anterior?
• ¿Qué parte de ella es solo casi siempre cierta?
• ¿Qué tiene todo eso que ver con GBM?

Está de moda usar caret para estas cosas, pero yo estoy todavía acostumbrado a hacerlas a mano. Creo, además, que es poco instructivo ocultar estas cuestiones detrás de funciones de tipo caja-negra-maravillosa a quienes se inician en el mundo de la construcción y comparación de modelos. Muestro, por tanto, código bastante simple para la validación cruzada de un modelo con R:

# genero ids
ids <- rep(1:10, length.out = nrow(cars))

# Nota: da igual si nrow(df) no es múltiplo de 10

# los aleatorizo
ids <- sample(ids)

# esto devuelve una lista de dfs:
preds.cv <- lapply(unique(ids), function(i){
  preds <- predict(lm(dist ~ speed,
    data = cars[ids != i,]), cars[ids == i,])
  data.frame(
    preds = preds,
    real = cars[ids == i,]$dist)
})

# "apilo" los dfs:
preds.cv <- do.call(rbind, preds.cv)

# calculo el rmse
rmse <- sqrt(mean((preds.cv$preds - preds.cv$real)^2))

Sí, estoy usando el RMSE aunque sea un detractor del mismo.

Doy clase en algunos máster de ciencia de datos. Estos máster suelen concluir con la realización de algún proyecto completo.

Ya sabemos cómo funcionan estas cosas en el medio académico: se busca cualquier cosa, se masomenos resuelve y se archiva. Sin recorrido ni impacto.

A mí me interesa proponer proyectos que tengan cierta trascendencia. El año pasado invité aun alumno a construir un sistema de predicción de plazas de aparcamiento disponibles en el sistema de bicicletas púbico de Zaragoza (dado que los datos están disponibles). El paquete MicroDatosEs mejoró como subproducto de otro proyecto. Etc.

La primera es la factorización positiva de matrices positivas. La otra, como bien titula la entrada, los procesos de Poisson autoexcitados.

Por eso no podía dejar de traer a la atención de mis lectores seismic. Aunque lo de Twitter ya huela.

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:

[N]NMF (se encuentra con una o dos enes) es una técnica de análisis no supervisado emergente. Se cuenta entre mis favoritas.

[N]NMF significa non negative matrix factorization y, como SVD, descompone una matriz M como UDV'. Solo que, en este caso, las entradas de M son todas positivas. Y la descomposición es UV', donde las entradas de ambas matrices son también positivas.

¿Qué tipo de matrices tienen entradas estrictamente positivas?

• Las resultantes de cuestionarios donde sujetos (filas) valoran (de 0 a 10) objetos, propuestas, etc. (columnas).
• Las que respresentan clientes (filas) que compran (un determinado número >= 0) de productos (columnas).
• …

Y acabo con un instrumento (el paquete NMF de R) y el análisis de una encuesta realizado con dicha técnica para que la veáis en acción.