Stan

Vivimos en un mundo opaco: como en los cuentecillos de Asimov, somos usuarios de tecnologías que ni conocemos ni controlamos. Parametrizamos nuestras máquinas y las echamos a correr. Poco más podemos hacer que fiarnos de quienes nos las proporcionan. Luego pasan cosas como que, de repente, resulta que Stan, en las últimas versiones, ha estado produciendo muestras sesgadas. ¿Qué resultados condicionará eso río abajo? Un caso mucho más famoso es el de la resonancia magnética (fMRI): un error en el software concomitante pone bajo sospecha hasta 40000 artículos sobre estudios del cerebro. Precisamente, por lo mismo. ...

Stan. Modelo multinivel. Variable categórica. Codificación con ceros y unos. Matriz. Coeficiente vector[n_ccaa] Cccaa. Sin priori. Catástrofe: (Coeficientes hasta 15000. Sin tasa, con tiempo. Los valores desorbitados, en ceros de la dummy). Priori. for (i in 1:n_ccaa) Cccaa[i] ~ cauchy(0, 20); ¿Por qué no? Tachán: (¿Para qué verbos?)

Por culpa de una nevera que no enfriaba como era debido, veinte años después, estoy repasando mi termodinámica: entropía, ciclo de Carnot, etc. Por culpa de Stan estoy repasando mi mecánica hamiltoniana. Y lo estoy disfrutando muchísimo. Dizque hay una exposición del Bosco en El Prado. Que si cuesta 16 euros. Que si solo puedes ver los cuadros de lejos porque hay toneladas de gente del extrarradio que hace su visita anual al centro. Pero, sobre todo, que es goce estético para los que no pueden apreciar otra cosa. Paso. Déjeseme hacer en paz lo que me gusta.

Continuando con la entrada anterior, ahora, números. Primero, el planteamiento (cuatro partidos, etc.): probs <- c(4, 3, 2, 1) probs <- probs / sum(probs) partidos <- letters[1:length(probs)] Nos hará falta más adelante library(plyr) library(rstan) library(ggplot2) library(reshape2) Sigo con el proceso de muestreo. Reitero: cada encuestador enseña al encuestado una tarjeta al azar donde aparece el nombre de dos partidos y le pregunta si ha votado (o piensa votar) a alguno de ellos. n <- 3000 resultados <- data.frame( tarjeta = sample(1:nrow(tarjetas), n, replace = T), partido = sample(partidos, n, prob = probs, replace = T)) resultados <- data.frame( tarjetas[resultados$tarjeta,], partido = resultados$partido) resultados$coincide <- resultados$partido == resultados$partido1 | resultados$partido == resultados$partido2 # proporciones reales en la muestra props.muestra <- table(resultados$partido) / nrow(resultados) # resultados agregados (por tarjeta) resultados.agg <- ddply( resultados, .(partido1, partido2), summarize, total = length(partido1), coincidencias = sum(coincide)) Y ...

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

Las diapositivas de mi charla Datos, modelos y parámetros en el grupo Machine Learning Spain pueden verse/bajarse de aquí.

Porque voy a dar una charla en él. Es este jueves, por la tarde, en el Campus de Google de Madrid (los detalles). Se tratará de una introducción a y justificación de aproximaciones más bayesianas de lo habitual a problemas reales del análisis de datos. Que comenzará con una explicación sobre cuándo 100% no significa 100% para terminar con lo que viene siéndome habitual últimamente: un ejemplo en rstan con su discusión.

y <- c(rnorm(1000), rnorm(2000, 1, 0.5)) es una mezcla de dos normales (N(0, 1) y N(1, 0.5)) con pesos 1/3 y 2/3 respectivamente. Pero, ¿cómo podríamos estimar los parámetros a partir de esos datos? Se puede usar, p.e., flexmix, que implementa eso del EM. Pero en el librillo de este maestrillo dice library(rstan) y <- c(rnorm(1000), rnorm(2000, 1, 0.5)) codigo <- " data { int<lower=1> K; // number of mixture components int<lower=1> N; // number of data points real y[N]; // observations } parameters { simplex[K] theta; // mixing proportions real mu[K]; // locations of mixture components real<lower=0> sigma[K]; // scales of mixture components } model { real ps[K]; // temp for log component densities sigma ~ cauchy(0,2.5); mu ~ normal(0,10); for (n in 1:N) { for (k in 1:K) { ps[k] <- log(theta[k]) + normal_log(y[n],mu[k], sigma[k]); } increment_log_prob(log_sum_exp(ps)); } }" fit <- stan(model_code = codigo, data = list(K = 2, N = length(y), y = y), iter=48000, warmup=2000, chains=1, thin=10) En el código anterior no sé si queda claro cómo cada punto $y_i$ sigue una distribución (condicionada a los parámetros) con densidad $\theta_1 \phi(y_i, \mu_1, \sigma_1) + \theta_2 \phi(y_i, \mu_2, \sigma_2)$. ...

Las diapositivas de mi charla sobre rstan en el grupo de usuarios de R de Madrid del 2016-02-11 están aquí. (Y los vídeos).

Este jueves (2016-02-11), a las 19:00, hablaré de rstan y de rstanarm en Medialab-Prado dentro de la reunión de usuarios de R de Madrid. Con el concurso de estos paquetes, replantearé tres problemas estadísticos conocidos desde una óptica bayesiana: Pruebas de hipótesis Regresión lineal Modelos estructurales de series temporales Si quieres asistir, reserva tu plaza aquí. Probablemente, discutiré todos esos modelos en estas páginas en los próximos días, además de colgar las diapositivas y sus fuentes.