Cuantiles, sí, pero ¿de qué tipo?

Porque resulta que los hay de varios tipos. En R, hasta nueve de ellos: set.seed(1234) muestra <- sort(rt(100, 3)) mis.cuantiles <- sapply(1:9, function(tipo) quantile(muestra, 0.834, type = tipo)) mis.cuantiles # 83.4% 83.4% 83.4% 83.4% 83.4% 83.4% 83.4% 83.4% 83.4% #0.9065024 0.9065024 0.8951710 0.8997036 0.9053693 0.9331290 0.9015846 0.9077920 0.9063154 Las definiciones de todos ellos pueden consultarse en Sample Quantiles in Statistical Packages. Las diferencias entre ellos, de todos modos, decrecen conforme aumenta el tamaño muestral: n.obs <- seq(100, 1e5, by = 1e3) res <- sapply(n.obs, function(n){ x <- rt(n, 3) diff(range(sapply(1:9, function(tipo) quantile(x, 0.834, type = tipo)))) }) plot(n.obs, log10(res), type = "l", xlab = "n obs", ylab = "discrepancia", main = "Diferencias entre los distintos tipos de cuantiles") ...

6 de marzo de 2017 · Carlos J. Gil Bellosta

Wikipedia + prophet

El otro día escribí sobre visitas a la Wikipedia. El otro día (posiblemente otro) oí hablar de prophet. Hoy con library(wikipediatrend) library(prophet) library(ggplot2) visitas <- wp_trend( "R_(lenguaje_de_programaci%C3%B3n)", from = "2010-01-01", to = Sys.Date(), lang = "es") mis.visitas <- visitas[, c("date", "count")] colnames(mis.visitas) <- c("ds", "y") pasado <- mis.visitas[1:1500,] m <- prophet(pasado) futuro <- make_future_dataframe(m, periods = nrow(mis.visitas) - 1500) prediccion <- predict(m, futuro) pred.plot <- plot(m, prediccion) pred.plot + geom_line(data = mis.visitas[1501:nrow(mis.visitas),], aes(x = ds, y = y), col = "red", alpha = 0.2) + xlab("fecha") + ylab("visitas") + ggtitle("Predicción de visitas a la página de R\nen la Wikipedia con prophet") construyo que muestra la predicción del número de visitas a la página de R en la Wikipedia basada en las primeras 1500 observaciones de la serie. Para las restantes, el gráfico muestra tanto los valores predichos como los reales (en rojo pálido). ...

3 de marzo de 2017 · Carlos J. Gil Bellosta

"Todas" las terrazas de Madrid

es un mapa en el que, en rojo, figuran todas (véase la coda) las terrazas de Madrid. Los datos están extraídos del censo de locales, sus actividades y terrazas de hostelería y restauración del ayuntamiento y están procesados con terrazas <- fread("http://datos.madrid.es/egob/catalogo/200085-17-censo-locales.txt") terrazas$coordenada_x_local <- as.numeric(gsub(",", ".", terrazas$coordenada_x_local)) terrazas$coordenada_y_local <- as.numeric(gsub(",", ".", terrazas$coordenada_y_local)) tmp <- terrazas[terrazas$coordenada_x_local > 1000, ] tmp <- terrazas[terrazas$coordenada_y_local > 3e6,] # UTM a siglo XXI library(rgdal) terrazas.utm <- SpatialPoints( cbind(tmp$coordenada_x_local, tmp$coordenada_y_local), proj4string=CRS("+proj=utm +zone=30")) terrazas.latlong <- spTransform(terrazas.utm, CRS("+proj=longlat")) library(ggmap) madrid <- get_map("madrid", zoom = 12) tmp <- as.data.frame(terrazas.latlong) ggmap(madrid) + geom_point( aes(x = coords.x1, y = coords.x2), data = tmp, size = .5, col = "red", alpha = 0.3) Sobre las cursivas de todas: ...

2 de marzo de 2017 · Carlos J. Gil Bellosta

Sobre una poco conocida y para nada menguante "brecha de género"

Con datos del INE sobre mortalidad he construido el gráfico que muestra las tasas de mortalidad relativas (la de hombres entre la de mujeres) desde 1975 para cada edad. Como no se aprecia debidamente el efecto que da pie a esta entrada, reorganizo los ejes (y promedio, ¡glups!, las tasas de mortalidad por grupos quinquenales de edad): Se observa una manifiesta tendencia creciente, uno de esos gender gaps, brechas de género o como quiera que se llamen a estas cosa en neolengua que, lejos de menguar, crece y crece. ...

1 de marzo de 2017 · Carlos J. Gil Bellosta

Consultando el número de visitas a páginas de la Wikipedia con R

Hace un tiempo probé el paquete wikipediatrend de R ya no recuerdo para qué. Desafortunadamente, el servicio que consulta debía de estar caído y no funcionó. Ahí quedó la cosa. Una reciente entrada de Antonio Chinchón en su blog me ha invitado a revisitar la cuestión y ahora, al parecer, stats.grok.se vuelve a estar levantado. Por lo que se pueden hacer cosas como: visitas <- wp_trend("R_(lenguaje_de_programaci%C3%B3n)", from = "2010-01-01", to = Sys.Date(), lang = "es") [Aquí ahorro al lector unos párrafos de pésima literatura.] ...

27 de febrero de 2017 · Carlos J. Gil Bellosta

Probando hunspell para el procesamiento de texto en español

El paquete hunspell de R permite procesar texto utilizando como soporte la infraestructura proporcionada por Hunspell, el corrector ortográfico que subyace a muchas aplicaciones en R. Existe una viñeta que ilustra el uso del paquete pero, como siempre, en inglés. En español las cosas son parecidas pero, como siempre, nunca exactamente iguales. En esta entrada, por lo tanto, voy a repasar partes de la viñeta aplicándolas a nuestra tan frecuentemente maltratada mas por ello no menos querida por algunos como yo (pausa) lengua. ...

20 de febrero de 2017 · Carlos J. Gil Bellosta

Una mala manera de perder un par de horas

Es esta: 156.67 * 100 # 15667 as.integer(156.67 * 100) #15666 Claro, hay que leer ?as.integer para enterarte de que, en realidad, la función que quieres usar es round. Una mala manera de perder un par de horas.

7 de febrero de 2017 · Carlos J. Gil Bellosta

El número efectivo de partidos

El número efectivo de partidos es el nombre de una página de la Wikipedia, que contiene la fórmula $$ N = \frac{1}{\sum_i p_i^2}$$ y excipiente alrededor. Aplicada a España (usando datos del CIS como proxy), Como casi siempre, el código: library(rvest) library(rvest) library(reshape2) library(plyr) library(zoo) url <- "http://www.cis.es/cis/export/sites/default/-Archivos/Indicadores/documentos_html/sB606050010.html" raw <- read_html(url) tmp <- html_nodes(raw, "table") tmp <- html_table(tmp[[2]], fill = TRUE) colnames(tmp)[1] <- "partido" tmp <- tmp[!is.na(tmp$partido),] tmp <- tmp[1:30,] tmp <- melt(tmp, id.vars = "partido") tmp <- tmp[tmp$value != ".",] tmp$value <- as.numeric(tmp$value) tmp$variable <- gsub("ene", "01-", tmp$variable) tmp$variable <- gsub("abr", "04-", tmp$variable) tmp$variable <- gsub("jul", "07-", tmp$variable) tmp$variable <- gsub("oct", "10-", tmp$variable) tmp$variable <- gsub("-0", "-200", tmp$variable) tmp$variable <- gsub("-1", "-201", tmp$variable) tmp$variable <- gsub("-9", "-199", tmp$variable) tmp$variable <- paste0("01-", tmp$variable) tmp$variable <- as.Date(tmp$variable, format = "%d-%m-%Y") dat <- tmp dat <- ddply(dat, .(variable), transform, total = value / sum(value)) res <- ddply(dat, .(variable), summarize, enp = 1 / (sum(total^2))) res <- zoo(res$enp, order.by = res$variable) plot(res, main = "Número efectivo de partidos\nen España(1996-2016)", xlab = "", ylab = "número efectivo de partidos")

25 de enero de 2017 · Carlos J. Gil Bellosta

Polinomios monótonos

Recibí un mensaje el otro día sobre polinomios monótonos. Mejor dicho, sobre el ajuste de datos usando polinomios monótonos. Frente a un modelo del tipo y ~ x (x e y reales) donde la relación entre las dos variables es manifiestamente no lineal y necesariamente monótina, p.e., creciente (por consideraciones previas), cabe considerar ajustar un polinomio monótono, i.e., realizar una regresión polinómica con la restricción adicional de que el polinomio de ajuste resultante sea monótono. ...

23 de enero de 2017 · Carlos J. Gil Bellosta

Va de si hay una o dos lambdas

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

18 de enero de 2017 · Carlos J. Gil Bellosta