R

Sutilezas de las licencias libres

R

Leyendo por ahí, he encontrado un comentario sobre el paquete RJSONIO de R en el que se recomendaba no usarlo por no ser libre.

El paquete, aparentemente, está liberado bajo una licencia BSD. Pero su pecado es que dentro de uno de los ficheros que contiene, src/JSON_parser.c, dice

The Software shall be used for Good, not Evil.

Más información, aquí.

No sé qué pensar sobre toda esta historia.

¿Quieres presentar algo en las Jornadas de Usuarios de R?

En varias de las ediciones de las Jornadas de Usuarios de R he formado parte del comité organizador, que se encarga, fundamentalmente, de la logística de la cosa. Este año, para variar, estoy en el comité científico.

Como integrante del cual, es labor mía tratar de animaros a que enviéis alguna propuesta de participación, que puede tener alguno de los siguientes formatos:

  • Una presentación de unos 20 minutos, mostrando alguna aplicación de R. Y no necesariamente en el mundo académico. Son también bienvenidas y apreciadas las aplicaciones en empresas e instituciones de  todo tipo. De hecho, una de las presentaciones más recordadas del año pasado, la de Antonio Sánchez Chinchón, trató de aplicaciones ludicomatemáticas de R.
  • Un taller (típicamente de 2 horas) para enseñar a utilizar alguna herramienta particular de R. En el pasado las ha habido de mapas, de gráficos, de Spark… ¡y recuerdo uno sobre plyr y reshape2 que impartí en 2010 cuando esos paquetes eran una novedad y una rareza!

Hay tiempo hasta el primero de mayo para realizar las propuestas. Los detalles pueden consultarse aquí.

Mezclas de distribuciones con Stan

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 $latex y_i$ sigue una distribución (condicionada a los parámetros) con densidad $latex \theta_1 \phi(y_i, \mu_1, \sigma_1) + \theta_2 \phi(y_i, \mu_2, \sigma_2)$.

Pequeño bug en ggmap: no pinta el último tramo de una ruta

R

Supongo que no debería escribirlo aquí sino comunicárselo a quien mantiene ggmap. Pero ya tuve una experiencia mejorable con él y dos no serán. Así que lo cuento por acá.

La mayor parte del mérito en el descubrimiento, en cualquier caso, es de una alumna de la clase de R que he dado hoy (en el momento en el que escribo, no en el que lees) en el Banco de Santander. No tengo su nombre ni tengo claro que quisiese que lo mencionase.

Ficheros KML con R y ggmap

Fácil:

library(maptools)
library(ggmap)

# un fichero bajado el Ayto. de Madrid
# (catálogo de datos abiertos)
rutas <- getKMLcoordinates("dat/130111_vias_ciclistas.kml")

# procesando el fichero kml
rutas <- lapply(1:length(rutas),
    function(x) data.frame(rutas[[x]], id = x))
rutas <- do.call(rbind, rutas)

# mapa de Madrid
mapa <- get_map("Madrid",
    source = "stamen", maptype = "toner",
    zoom = 12)

# pintando los tramos sobre el mapa
ggmap(mapa) + geom_path(aes(x = X1, y = X2,
    group = id), data = rutas,
    colour = "red")

produce

rutas_ciclistas_madrid

Nota: KML es esto.

Validación cruzada en R

Está de moda usar caret para estas cosas, pero yo estoy todavía acostumbrado a hacerlas a mano. Creo, además, que es poco instructivo ocultar estas cuestiones detrás de funciones de tipo caja-negra-maravillosa a quienes se inician en el mundo de la construcción y comparación de modelos. Muestro, por tanto, código bastante simple para la validación cruzada de un modelo con R:

# genero ids
ids <- rep(1:10, length.out = nrow(cars))

# Nota: da igual si nrow(df) no es múltiplo de 10

# los aleatorizo
ids <- sample(ids)

# esto devuelve una lista de dfs:
preds.cv <- lapply(unique(ids), function(i){
  preds <- predict(lm(dist ~ speed,
    data = cars[ids != i,]), cars[ids == i,])
  data.frame(
    preds = preds,
    real = cars[ids == i,]$dist)
})

# "apilo" los dfs:
preds.cv <- do.call(rbind, preds.cv)

# calculo el rmse
rmse <- sqrt(mean((preds.cv$preds - preds.cv$real)^2))

Sí, estoy usando el RMSE aunque sea un detractor del mismo.

Hoy se ha anunciado la propuesta de nueva página de la Comunidad R Hispano

R

Acabo de escribir a los socios de la Comunidad R Hispano acerca de la existencia de una propuesta para renovar la página de la asociación. Podéis ver la versión actual y la propuesta.

Logo de la Comunidad R Hispano

(Y agradezco muchísimo el trabajo de Paula López Casado, responsable de que la nueva página tenga un aspecto infinitamente más atractivo que la que lees).

(Además, esta entrada incluye el que será el nuevo logo de la Comunidad R Hispano).

storr: como Redis, pero con R

R

Probablemente no habéis utilizado nunca Redis. Redis es un sistema de almacenamiento basado en parejas clave-valor. Es similar a un diccionario de Python o a un entorno en R. Salvo que el almacenamiento es externo al proceso: los datos se guardan en un sistema distribuido y potencialmente ilimitado en cuanto a capacidad.

Si queréis probar algo parecido, además de los diccionarios y los entornos, podéis probar con storr , un paquete reciente de R. Aquí tenéis una minisesión de ejemplo:

Premoniciones de Tirole sobre sobre el R Consortium

R

A J. Tirole tiene Nobel de economía. En 2002 escribió un artículo, Some Simple Economics of Open Source, en el que trataba de explicar desde un punto de vista económico y de organización industrial el porqué de esa rareza. Aparte de cuestiones como si sería extrapolable a otros sectores distintos del del desarrollo de software.

En la sección sobre la reacción de las compañías de software frente al fenómeno del software libre tiene un apartado titulado viviendo simbióticamente de [no con] un proyecto de código abierto que termina con la frase (mi traducción):

Comparaciones de tres grupos: pruebas vs modelos

Una pregunta reciente en r-help-es se refería a la comparación en R de las proporciones en tres grupos. Obviando algunas pequeñas complicaciones en el problema, la respuesta canónica podría ser esta:

total <- c(56, 49,51)
positivos <- c(14, 10, 17)
prop.test(tmp$positivos, tmp$positivos + tmp$negativos)

# 3-sample test for equality of proportions without continuity correction
#
# data:  tmp$positivos out of tmp$positivos + tmp$negativos
# X-squared = 2.2289, df = 2, p-value = 0.3281
# alternative hypothesis: two.sided
# sample estimates:
#   prop 1    prop 2    prop 3
# 0.2500000 0.2040816 0.3333333

Los grupos no parecen ser desiguales.

Análisis estadístico de respuestas ocultas en encuestas

A veces se hacen encuestas sobre temas sobre los que los encuestados son reticentes a revelar la verdad (p.e., ¿es Vd. un zombi?). Un procedimiento conocido para recabar tal tipo de información es el siguiente:

  • Se le invita al encuestado a tirar al aire una moneda con las caras etiquetadas con y no; la moneda no es una moneda porque tiene una probabidad conocida (y distinta del 50%) de caer en .
  • El encuestado responde sí si la respuesta a la pregunta y el resultado de la tirada de la moneda coinciden y no en caso contrario.

A partir de la proporción de respuestas positivas y conocida la probabilidad del de la moneda, $latex q$, es posible estimar la proporción $latex \theta$ de respuestas positivas a la pregunta de subyacente de interés en la muestra. Efectivamente, los síes tienen una distribución binomial $latex B(p) = B(q\theta + (1-q)(1-\theta))$ y, una vez estimado (por máxima verosimilitud) $latex \hat{p}$, puede despejarse $latex \hat{p}$ de $latex \hat{p} = q\hat{\theta} + (1-q)(1-\hat{\theta})$ para obtener