Hace un tiempo resumí los GBMs (Gradient Boosting Machines) en una línea. Hoy comienzo una serie de varias entradas para que nadie tenga excusa de no saber de qué va la cosa. Arranco con una mentira sugerente. Porque lo que voy a contar no es del todo cierto, pero motiva lo que vendrá después.

Consideremos un conjunto de datos medio famoso: el de los precios de los alquileres en Múchich. Comencemos con un modelo sencillo, una regresión lineal que relacione el precio del alquiler con los metros cuadrados, i.e.,

O tal dice lo que expongo a continuación.

Paquetes necesarios:

library(rvest)
library(caRtociudad)
library(reshape2)
library(ggmap)
library(plyr)
library(TSP)

Extracción de las provincias y sus capitales (de la Wikipedia):

capitales <- read_html("https://es.wikipedia.org/wiki/Anexo:Capitales_de_provincia_de_Espa%C3%B1a_por_poblaci%C3%B3n")
capitales <- html_nodes(capitales, "table")
capitales <- html_table(capitales[[1]])$Ciudad

capitales <- capitales[!capitales %in%
  c("Las Palmas de Gran Canaria",
    "Melilla", "Ceuta", "Mérida",
    "Santa Cruz de Tenerife",
    "Santiago de Compostela",
    "Palma de Mallorca")]

Y sus coordenadas:

coordenadas <- ldply(capitales, function(x) {
    tmp <- cartociudad_geocode(x)[1,]
    res <- data.frame(ciudad = x, provincia = tmp$province, lat = tmp$latitude, lon = tmp$longitude)
    if(is.na(res$lat)){
      tmp <- geocode(paste(x, "España"))
      res$lat <- tmp$lat
      res$lon <- tmp$lon
    }
    res
  })

# Pobre Logroño: ¡Cartociudad lo ubica en Asturias!
coords.logrono <- geocode("Logroño")
coordenadas$lat[coordenadas$ciudad == "Logroño"] <- coords.logrono$lat
coordenadas$lon[coordenadas$ciudad == "Logroño"] <- coords.logrono$lon

Construcción de la matriz simétrica de distancias (¡tarda un buen rato!):

Una de las cosas de las que me acuerdo de cuando leía es un parrafito de Mi idolatrado hijo Sisí en el que Delibes ponía en la boca no sé si de alguno de sus protagonistas o del narrador en el que se daba cuenta de la anormalidad histórica que supuso el tiempo de la II República: de repente, la gente hacía cosas que nunca había hecho y que nunca había vuelto a hacer: hablar a todas horas y con todo el mundo de política. La política entraba en los círculos de amigos, en la sobremesa de las familias, etc.

Tengo la sensación de que un lenguaje funcional (como Scala) está particularmente bien adaptado al tipo de operaciones que exige MCMC.

Juzguen Vds.

Primero, genero datos en R:

datos <- rnorm(500, 0.7, 1)
writeLines(as.character(datos), "/tmp/datos.txt")

Son de una normal con media 0.7. En el modelo que vamos a crear, suponemos conocida (e igual a 1) la varianza de la normal y trataremos de estimar la media suponiéndole una distribución a priori normal estándar. Y con Scala, así:

Aunque esta entrada es sin duda resabida de los más de mis lectores, quedarán los que aún no sepan que ciertas distribuciones no tienen media. Condición necesaria para que una distribución la tenga es que

$$ \int_{-\infty}^\infty |x| f(x) dx$$

tenga un valor finito, cosa que, por ejemplo, no cumple la de Cauchy. Igual hay a quien esto le parece una rareza matemática, un entretenimiento de math kiddies sin implicaciones prácticas. Además, porque para que que la integral anterior diverja se necesita que las distribuciones puedan tomar valores arbitrariamente altos y las que se manejan en la práctica están acotadas si no por el número de átomos del universo por el de céntimos de bolívar venezolano necesarios para comprar todas las cosas que caben en el ancho mundo.

Hay platos con nombre. P.e., tortilla de patata o tiramisú. También hay distribuciones (de probabilidad) con nombre. P.e., normal, binomial, Poisson, hipergeométrica.

Hay quienes quieren saber (1) todas (o muchas) de esas distribuciones con nombre y (2), dados unos datos, cuál de ellas siguen. Esta entrada va a tener la url a la que de ahora en adelante remita a quien me las formule.

A pesar de que algunos platos tienen nombre, el otro día se podía probar en el Diverxo espárrago blanco a la mantequilla negra con emulsión de leche de oveja, espardeña y salmonete. Que no es ni tortilla de patata, ni tiramisú ni otra cosa con nombre que se le parezca.

(Aviso: esta entrada podría competir dignamente en una competición de titulares engañosos. Es posible que si no sepas de qué hablo regularmente te interese más esto).

En España hay pruebas de acceso a la universidad que y en algunos sitios publican las notas de corte para acceder a determinados estudios. Las he bajado escrapeando El País así

library(rvest)
library(plyr)
library(rstan)
library(reshape2)

options(mc.cores = 2)

url <- "http://elpais.com/especiales/universidades/"

pagina     <- read_html(url, encoding = "UTF8")
urls_provs <- html_nodes(pagina, "a")
urls_provs <- html_attr(urls_provs, "href")
urls_provs <- paste0("http://elpais.com", urls_provs[grep("centro/provincia", urls_provs)])

foo <- function(url){
  tmp  <- read_html(url)
  urls <- html_nodes(tmp, "a")
  urls <- html_attr(urls, "href")
  paste0("http://elpais.com", urls[grep("^/especiales/universidades/titulacion/universidad", urls)])
}

urls_univs <- sapply(urls_provs, foo)
urls_univs <- unique(unlist(urls_univs))


foo <- function(url){
  tmp <- read_html(url)
  lugares <- html_nodes(tmp, xpath = "//*/div[@class = 'lugar']")
  data.frame(carrera = html_text(html_nodes(lugares, xpath = "//*/a[@class = 'carrera']")),
              sede = html_text(html_nodes(lugares, xpath = "//*/p[@class = 'escuela']/span/a")),
              nota = html_text(html_nodes(tmp, xpath = "//*/div[@class='nota']/p/text()")))
}

res <- ldply(urls_univs, foo)

notas <- res
notas$nota <- as.numeric(as.character(res$nota))

# limpieza de datos
notas$nota[notas$nota > 5000] <- notas$nota[notas$nota > 5000] / 1000
notas <- notas[notas$nota > 0,]

notas <- notas[order(notas$nota), ]

con el objetivo de estudiar el efecto de la universidad / sede y de la carrera en el punto de corte. Esencialmente, quiero hacer algo así como lmer(nota ~ 1 + (1 | sede) + (1 | carrera), data = notas), pero hay una complicación: como creo que mis lectores sabrán, las notas de acceso tienen un valor mínimo, el del aprobado, 5. Eso significa que, de alguna manera, están censuradas por la izquierda. El modelo resultante es algo así como

Porque ahí puedes tomarte una foto tal que

o

y luego tuitear cosas como

After #StrataHadoop: @jayusor and me in front of Bayes grave with my new @OReillyMedia book (emulating @gilbellosta) pic.twitter.com/SlguJIeLw0

— Antonio Sánchez Chinchón (@aschinchon) June 4, 2016

Para mayor referencia (y por tenerlo a mano cuando vuelva),

Me estoy volviendo intolerante al ruido. Y esta mañana (¿qué carajos hago levantado tan temprano?) no había forma de que dejase de sonar la alarma de unos andamios de la plaza, no paraba la batidora del bar desde donde escribo y, encima, esto, esto, esto, esto, esto, esto,…

Son todas noticias relacionadas con la publicación de esto, un artículo que describe un estudio clínico (¡con 84 sujetos!) en el que se comparan dos grupos (uno tratado y otro no) que,

La versión larga, en inglés y capada para quienes tengáis bloquador de publicidad es esta. La versión abreviada para el hombre ocupado de hoy en día (y con mis comentarios) es:

• Un tipo con mucha vocación estudia física y acaba en el CERN
• Trabaja en lo de la partícula de Higgs y obtiene peor de los resultados posibles: la partícula existe tal cual predecía la teoría
• Lo de la particulita queda ahí porque, desafortunadamente, no hay nada más que rascar
• Se acaban las subvenciones, las becas, etc. y todos los físicos vocacionales tienen que buscarse la vida
• El protagonista de nuestra historia fabrica apps de telefonillo para carreras de caballos

Mis comentarios: