R

Un pequeño problema de probabilidad

El tuit

de John Allen Paulos me indujo a escribir

number.numbers <- function(n){
  sum(cumsum(sample(0:n)) < n) + 1
}

res <- replicate(10000, number.numbers(1000))

código con el que, efectivamente, puede comprobarse que la media es, efectivamente, e.

Ahora bien, ¿alguien se atreve a explicar por qué?

(No leas esta pista: (s??)?s??).

rPython, ya en Windows

R

Aprovechando que por un lado las circunstancias han querido que ahora disfrute de más tiempo libre; que, por otro, mi paquete rPython parece ir ganando aceptación y, finalmente, que tengo varios correos pendientes clamando por una versión en Windows, he pasado unos ratos tratando de hacer el proceso de instalación lo menos pesado y manual que me ha sido posible.

Y el resultado ha sido este.

Así que si alguien todavía sigue usando Windows y tiene interés en interactuar con Python desde R (aunque solo sea por aburrimiento), que lo pruebe. Y si algo se rompe, que me lo haga saber para parchear el proceso.

La red Asia

La red Asia es esto:

Es decir, una red bayesiana. Una red bayesiana clásica sobre la que los interesados podrán saber más leyendo lo que Lauritzen y Spiegelhalter dejaron escrito sobre ella en 1988.

Pero la idea básica es la siguiente:

  • Los nodos superiores (visita a Asia, fumador) son variables observables sobre el comportamiento de unos pacientes.
  • Los nodos inferiores (rayos X, disnea) son variables también observables, síntomas de esos pacientes.
  • Los nodos centrales, los más importantes, no son observables: son diversas enfermedades que pudieran estar padeciendo los individuos en cuestión.

La pregunta que ayuda a resolver esta red bayesiana es la siguiente: conocidas (¡o no!) las variables observadas, ¿cuál es la probabilidad de que un paciente dado padezca alguna de las enfermedades (tuberculosis, bronquitis o cáncer de pulmón) correspondientes a los nodos centrales?

5000 paquetes de R en CRAN

R

CRAN (el Comprehensive R Archive Network) acaba de superar la barrera de los 5000 paquetes. Reproduzco (y traduzco) un mensaje de Henrik Bengtsson en la lista de desarrollo de R que ofrece más información al respecto:

Pasar de los 4000 a los 5000 paquetes tardó 14.5 meses, un paquete nuevo cada 10.5 horas. Detrás de cada paquete hay gente real. En total, 2900 personas mantienen paquetes, de los que 350 son nuevos, y muchos más que contribuyen a ellos. Mucha más gente a la que agradecer que nunca antes (¡y a quienes hay que citar adecuadamente en las publicaciones!).

Importancia de variables en árboles

Los árboles (o árboles de inferencia condicional) valen fundamentalmente para hacerse una idea de cómo y en qué grado opera una variable en un modelo controlando por el efecto del resto. Su valor reside fundamentalmente en la interpretabilidad.

No obstante lo cual, no es infrecuente construir árboles muy grandes. Y el tamaño dificulta censar qué variables y en qué manera aparecen. Por eso me vi obligado recientemente a crear un pequeño prototipo para extraer el peso de las variables de un árbol.

Un récord personal

El otro día, casi por error, cargué este dataframe en R:

dim(raw)
# [1] 115318140         4

Es todo un récord personal logrado en un servidor con 24GB de RAM bastante caro.

El anterior estaba en otro de algo así como 20 millones de filas y unas 6 o siete columnas. Eso sí, logrado en tiramisu, mi ordenador personal de 8GB de RAM de 400 euros (monitor incluido).

Os preguntaréis si pude hacer algo con ese monstruo. La verdad es que sí: pude muestrear un 10% de las filas y trabajar con ellas sin problemas.

Contando hexágonos en paralelo

Dicen que para realizar gráficos de dispersión con muchos datos no es desaconsejable usar celosías hexagonales. Por motivos que no vienen al caso, me interesa poder realizarlas en paralelo.

El código disponible en R (hexBinning de fMultivar o el de geom_hex de ggplot2) es feo, ininteligible y, en particular, no es paralelizable (o mapreducible). No lo es porque cada hilo, por diseño del algoritmo, crea hexágonos excéntricos.

Así que he desarrollado un algoritmo para crear celosías hexagonales paralelizable. Además, creo que es algo más inteligible que los dos mencionados y, me temo, igual de feo. Pero vectorizado, eso sí (es decir, sin un maldito bucle). Es así:

Sexo, deporte y la cantidad de información mutua

Perdón por el titular. No soy inasequible a las modas.

La cuestión del día de hoy es la siguiente: tenemos una variable X inobservable y otra variable Y potencialmente correlacionada con X. ¿Cuánto podemos decir de X de conocida Y?

Supongamos que ambas son binarias. Si conozco Y poseo 1 bit de información. Si solo conozco X (que me da pistas sobre Y) conoceré una fracción de un bit de información (sobre Y).

Recordatorio: V Jornadas de Usuarios de R, diciembre 2013, Zaragoza

Que sirva esta entrada de recordatorio de que la comunidad de usuarios de R va a celebrar sus V Jornadas en Zaragoza este mes de diciembre.

Este año, como aliciente adicional, contamos con no uno sino tres concursos con el patrocinio de Synergic Partners y Telefónica Digital.

Como miembro del comité organizador de las jornadas, reitero mi invitación a quienes siguen estas páginas a que acudan, participen y, en el peor de los casos, ayuden a difundirlas entre quienes consideren interesados en ellas.

El cuarteto de Anscombe

F. Anscombe escribió en 1973 el artículo Graphs in Statistical Analysis para subrayar la importancia de los gráficos en el análisis estadístico.

Esencialmente, el artículo se limita a presentar cuatro conjuntos de datos distintos con la misma media, varianza, correlación y recta de regresión (excúsenseme los abusos del lenguaje). Sin embargo tienen aspectos muy distintos:

¿Os interesa saber más al respecto? Pues tenéis esto para satisfacer el apetito de culturilla y esto otro para jugar con el cuarteto en R.

Medianas ponderadas en R

La mediana de 1:3 es 2. Pero puede ser que queramos dar a 1:3 los pesos 2, 1, 2. En ese caso, el cálculo de la mediana sigue siendo sencillo (y sigue siendo 2). Pero la situación puede complicarse más.

Mientras los pesos sean enteros, todavía pueden usarse trucos:

x <- 1:3
pesos <- c(2,1,2)
median(rep(x, times = pesos ))

¿Pero qué hacemos cuando hay pesos fraccionarios? Bueno, en realidad, podemos ordenar:

n <- 1000

x <- runif(n)
pesos <- runif(n)
o <- order(x)
x.o <- x[o]
pesos.o <- pesos[o]
x.o[min(which(cumsum(pesos.o) > .5 * sum(pesos.o)))]

Pero me parece más limpio usar el paquete quantreg:

Las V Jornadas de Usuarios de R, en Zaragoza

Escribo para anunciar públicamente que están en marcha las V Jornadas de Usuarios de R. Se celebrarán este año en Zaragoza, los días 12 y 13 de diciembre.

Todavía no está disponible el programa (que, en cierto modo, es responsabilidad de vosotros: estáis invitados a enviar propuestas de ponencias y talleres). Tenemos un concurso cuyas bases podrían todavía modificarse si un generoso patrocinador asumiese su financiación.

Y eso, que estáis todos invitados a esta nueva edición de las jornadas.

¿Nos ayudáis a mejorar r-es.org?

R

Hora es cumplida, creo yo, de repensar el portal de la Comunidad R Hispano. Así que he pensado en pulsar el criterio —que estimo sobremanera— de mis visitantes y solicitar de ellos (vía comentario a esta entrada) sugerencias. Por acotar el tema, sugiero que vayan encaminadas a dirimir estas dos cuestiones:

  1. ¿Cuál debería ser el objetivo de un portal de esas características?
  2. ¿Cómo debería organizarse para alcanzar mejor esos objetivos?

Y como colofón, ¿conocéis algún modelo aplicable y que funcione?