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

Queríamos desentrañar los misterios de las partículas subatómicas y obtuvimos una app para las carreras de caballos

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

7 de junio de 2016 · Carlos J. Gil Bellosta

Acceso a Google Analytics desde R

Google Analytics puede usarse desde su consola o bien descargando datos y procesándolos por tu cuenta. Para lo cual, desde R, require(RGoogleAnalytics) client.id <- "1415926535-u377en6un7lugar2de7lamancha0de1cuyo5nombre0m.apps.googleusercontent.com" client.secret <- "CEcI5nEst6pAs6Un2SecREt6-f8nt" token <- Auth(client.id,client.secret) #save(token,file="~/.ga_token_file") Obviamente, para lo anterior: Hay que instalar y cargar los paquetes relevantes Tienes que usar tu propio id y secreto de cliente como indica aquí Tienes que tener una cuenta en Google Analytics, claro Además, puedes descomentar la última línea si quieres guardar tus credenciales para futuros usos (con las debidas medidas de seguridad). Tras lo cual, ...

6 de junio de 2016 · Carlos J. Gil Bellosta

R sobre el EC2 de Amazon hace casi siete años: una concesión a la melancolía

Corría el año 2009 cuando comencé mi segunda aventura bloguera (nadie, yo incluido, quiere rememorar la primera) cuando Raúl Vaquerizo tuvo la caridad de aceptarme como colaborador en Análisis y Decisión. En diciembre de aquel año escribí cómo utilizar R en una cosa que entonces comenzaba a sonar: la nube y, en concreto, el servicio EC2 de Amazon. El resultado, probablemente totalmente desfasado, fue este. Material de hemeroteca, alimento de melancolías.

3 de junio de 2016 · Carlos J. Gil Bellosta

Detección de "outliers" locales

Aunque outlier local parezca oxímoron, es un concepto que tiene sentido. Un outlier es un punto dentro de un conjunto de datos tan alejado del resto que diríase generado por un mecanismo distinto que el resto. Por ejemplo, puedes tener las alturas de la gente y alguna observación que parece producto de otra cosa como, por ejemplo, errores mecanográficos en la transcripción. Un outlier está lejos del resto. Pero, ¿cuánto? Con ciertas distribuciones tiene sentido pensar que los outliers son puntos a una distancia superior a nosecuántas desviaciones típicas de la media. Más en general, fuera de un determinado círculo. Una medida similar: serían outliers aquellos puntos que a una determinada distancia solo tienen un determinado porcentaje (pequeño) del resto. Todas estas son medidas globales. ...

2 de junio de 2016 · Carlos J. Gil Bellosta

¿Alguien podría identificar tirios y troyanos?

Con los datos pcts <- cbind( c(35.7, 19.6, 6.6, 16.6, 9.6), c(0.3, 0.2, 0.2, 0.3, 0.8), c(25.0, 14.9, 10.7, 32.7, 12.9), c(1.6, 8.0, 8.5, 6.5, 7.9), c(11.0, 18.7, 7.9, 12.7, 8.0), c(3.2, 21.5, 52.9, 16.7, 47.9) ) totales <- c(1102, 975, 596, 638, 174) tabla <- round(t(pcts * totales / 100)) y el concurso de library(MASS) biplot(corresp(tabla, nf = 2)) genero que a lo mejor no resulta demasiado interesante si no añado que las columnas se refieren a partidos políticos y las filas a cadenas en las que, según el CIS, sus votantes prefieren para seguir la actualidad política. Eso sabido, ¿cuál es cuál?

1 de junio de 2016 · Carlos J. Gil Bellosta

El extraño caso de la media empírica menguante

La distribución lognormal es la exponencial de una distribución normal. Su media, Wikipedia dixit, es $\exp(\mu + \sigma^2 /2)$. Dada una muestra de la distribución lognormal (y supuesto, por simplificar, $\mu=0$), podemos calcular su media y una estimación de su $\sigma$ y calcular $\exp(\sigma^2 /2)$ y uno pensaría que los valores deberían ser similares. Mas pero sin embargo, library(ggplot2) set.seed(123) sigmas <- seq(1, 10, by = 0.1) res <- sapply(sigmas, function(sigma){ a <- exp(rnorm(1e6, 0, sigma)) mean(a) / exp(var(log(a))/2) }) tmp <- data.frame(sigmas = sigmas, medias = res) ggplot(tmp, aes(x = sigmas, y = medias)) + geom_point() + geom_smooth() produce ...

31 de mayo de 2016 · Carlos J. Gil Bellosta