GBM (II): Minimización de funciones, pérdidas cuadráticas, residuos y gradientes

Para minimizar una función $\phi(x)$ es habitual utilizar un procedimiento iterativo: a partir de un punto inicial $x_0$ se salta a $x_1 = x_0 - \lambda \nabla \phi(x_0)$ (donde $\lambda$ es un número pequeño predefinido de antemano), y de ahí, sucesivamente, a $$ x_n = x_{n-1} - \lambda \nabla \phi(x_{n-1}).$$ Porque, típicamente, como cuando uno está en el monte y da un paso corto en la dirección opuesta a la de máxima pendiente, sucede que ...

22 de junio de 2016 · Carlos J. Gil Bellosta

GBM (I): Una mentira sugerente

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únich. Comencemos con un modelo sencillo, una regresión lineal que relacione el precio del alquiler con los metros cuadrados, i.e., ...

21 de junio de 2016 · Carlos J. Gil Bellosta

6602.767 km alrededor de España para visitar todas sus capitales de provincia

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

20 de junio de 2016 · Carlos J. Gil Bellosta

Evolución histórica de la deuda del ayuntamiento de Madrid

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

17 de junio de 2016 · Carlos J. Gil Bellosta

Metropolis-Hastings en Scala

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

16 de junio de 2016 · Carlos J. Gil Bellosta

Distribuciones sin media: ¿qué pueden suponer en la práctica?

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

15 de junio de 2016 · Carlos J. Gil Bellosta

Las distribuciones (y platos) con nombre

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

14 de junio de 2016 · Carlos J. Gil Bellosta

Censura a la izquierda en las universidades españolas

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

13 de junio de 2016 · Carlos J. Gil Bellosta

Si vas a Londres, déjate caer por (51.523841, -0.089310)

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),

10 de junio de 2016 · Carlos J. Gil Bellosta

Ruido de alarmas, ruido de p-valores; mucho, mucho ruido, tanto, tanto ruido

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

8 de junio de 2016 · Carlos J. Gil Bellosta