Stan

Envalentonado por el comentario de Iñaki Úcar a mi entrada del otro día, que me remitía a este artículo, decidí rizar el rizo y crear intervalos de confianza no ya discontinuos sino con otra propiedad topológica imposible: homeomorfos con un toro. Y aquí está: El modelo, el código y demás, library(rstan) library(ggplot2) n <- 100 a1 <- 1 a2 <- 1 sigma <- 0.4 datos <- data.frame(x1 = rnorm(n, 2, 0.1), x2 = rnorm(n, 2, 0.1)) datos$y <- a1^datos$x1 + a2^datos$x2 + rnorm(n, 0, sigma) codigo <- " data { int<lower=1> N; real y[N]; real x1[N]; real x2[N]; } parameters { real<lower=-3, upper="3"> a1; real<lower=-3, upper="3"> a2; real<lower=0, upper="3"> sigma; } model { for (n in 1:N) y[n] ~ normal(fabs(a1)^x1[n] + fabs(a2)^x2[n], sigma); }" fit <- stan(model_code = codigo, data = list(N = length(datos$y), y = datos$y, x1 = datos$x1, x2 = datos$x2), iter=40000, warmup=2000, chains=1, thin=10) res <- as.data.frame(fit) ggplot(res, aes(x = a1, y = a2)) + geom_point(alpha = 0.1) De nuevo, no son intervalos propiamente dichos, lo convengo. Pero son configuraciones más fieles al espíritu de lo que un intervalo de confianza es y representa que su(s) letra(s) I N T E R V A L O. ...

Un año, el 2016, mueren 1160 personas en accidentes de tráfico. El anterior, 1131, i.e., 29 menos. Ruido estadístico aparte, ¿aumentan? Comenzamos a optar. Primera elección subjetiva: son muestras de una Poisson de parámetro desconocido. La pregunta: ¿el mismo? Una manera de estudiar lo anterior es plantear 1160 ~ poisson(lambda * (1 + incr)) 1131 ~ poisson(lambda) y estudiar la distribución de incr. Que a saber qué distribución tendrá (teóricamente). Pero, ¿importa? Mejor que rebuscar a ver qué distribución podría tener la cosa, basta con envolverlo en un poco de seudo-C++, ...

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