R

Hoy comienzo a enseñar un curso de introducción a la programación para recién graduados que comenzarán un máster de matemáticas aplicadas con incursiones en la llamada ciencia de datos. Serán 4 sesiones con el siguiente contenido: Sesión 1, programación imperativa: variables, condicionales y bucles. Sesión 2, programación orientada a objetos. Sesión 3, colecciones: listas, tuplas, conjuntos, diccionarios, etc. Sesión 4, programación funcional: map, reduce, fold, foldLeft, scan, filter, etc. Los lenguajes a utilizar serán R y Python (via Jupyter). No me he atrevido a añadir Scala (como ejemplo de cómo deben hacerse las cosas, además de ser un lenguaje, para variar, tipado y no interpretado) por falta de tiempo. ...

Desde el principio de mis tiempos he seleccionado variables relevantes como subproducto de los árboles primero y de los bosques aleatorios después. Cierto que he hecho casi inconfesables incursiones en los métodos stepwise, pero han sido marginales y anecdóticas. La idea es casi siempre la misma, se haga a mano o con ayuda de paquetes ad hoc: las variables importantes tienden a aparecer en el modelo (o submodelos), las otras no. Todo se reduce a contar y ponderar. Hay que discurrir un poco más cuando se sospecha (o consta) que existen variables altamente correlacionadas. ...

En los últimos tiempos se ha puesto de moda un subgénero periodístico que es una manera de generar artículos de acuerdo con el siguiente algoritmo: Se toma una lista de personas. Se cuenta en ella el número de mujeres (a) y de hombres (b). Si a >= b, GOTO 1; si no, se copipega y se mutatismutandea un manido argumento. No sabiéndome sustraer al encanto del último grito, he escrito y corrido ...

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: ...

A esto:

Siempre (aténganse los puristas al contexto) recomiendo comenzar con un árbol de decisión para, sobre esa base, ensayar métodos más potentes. Sobre todo si la precisión conviene más que la interpretabilidad. En la dirección opuesta se sitúan los árboles rápidos y frugales. Un árbol rápido y frugal es un tipo de árbol de decisión tal como La restricción que satisface (a diferencia de los árboles de decisión más habituales) es que: ...

Los productos de Apple, aun admitiendo su calidad, resuelven problemas que yo hace años que no tenía. Tanto data.table como dplyr vinieron a resolver problemas a los que muchos nos enfrentábamos con sudor y lágrimas. Ha aparecido recientemente una serie de tutoriales sobre ambos paquetes que recomiendo: El de data.table El de dplyr (parte I, parte II) Y mis comentarios: Para el 99% de mis problemas de manipulación de datos, me sobra con, además de R base, reshape2 y plyr. Para datos más grandes, me decanto por data.table. En gran medida, porque es previo a dplyr. No obstante, tengo la sensación de que dplyr acabará llevándose el gato al agua: tengo suficientes años como para haber presenciado sin que me las cuenten batallas anteriores: Beta vs VHS, WordPerfect vs Word, etc.

Jan Gorecki ha resumido las soluciones a las cincuenta preguntas más populares sobre el paquete data.table de R en Stack Overflow y las ha resumido en forma de tutorial aquí. Muy recomendable. Muy recomendable también data.table. Aunque me temo que el hadleyverse, y por razones que nada tienen que ver con la calidad de la cosa, no van a dejar de él, a medio plazo, ni las raspas.

rio es otro de esos desasosegantes paquetes de R. rio contiene esencialmente tres funciones, import, que lo lee todo export, que lo escribe todo y convert, que transforma un fichero de un formato a otro. Según su documentación, uno puede hacer cosas como export(mtcars, "mtcars.csv") export(mtcars, "mtcars.rds") export(mtcars, "mtcars.sav") para guardar mtcars en cualquiera de los formatos indicados por la extensión y luego x <- import("mtcars.csv") y <- import("mtcars.rds") z <- import("mtcars.sav") para importarlos sin complicaciones de separadores, encabezamientos, etc. Además, ya puestos, ...

Esto que cuento hoy puede ser muy útil: cómo mejorar los GLMs mediante la introducción de prioris (casi) a voluntad sobre los coeficientes. Usando el paquete arm de R, claro. De momento y porque aún tengo sucios los datos sobre los que me gustaría aplicar el modelo, extraeré un ejemplo de la ayuda de la función principal del paquete, bayesglm. Primero, preparo unos datos: n <- 100 x1 <- rnorm (n) x2 <- rbinom (n, 1, .5) b0 <- 1 b1 <- 1.5 b2 <- 2 y <- rbinom (n, 1, invlogit(b0+b1*x1+b2*x2)) Comenzamos con un glm de toda la vida. ...