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

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.

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

Nota: KML es esto.

Hago pública por su interés (parte de) una respuesta de Ramón Díaz Uriarte a un correo mío en el que yo sugería

que una vez que sabes especificar un modelo probabilístico para unos datos, p.e.,

• para la regresión lineal, y ~ N(a0 + a1 x1 +..., sigma)),
• para el test de Student, y0 ~ N(mu, sigma); y1 ~ N(mu + delta, sigma),
• etc. no hace falta saber qué es lm, ni el test de Student, ni nada. Cero teoría; sobre todo, de teoría tipo recetario. Se especifica el modelo (con una determinada sintaxis), se deja correr la cosa y a interpretar.

Su respuesta:

Tengo un grafo, g cuyas aristas pueden ser cualquier cosa susceptible de contaminarse. Me pregunto si la contaminación puede contagiarse a través del grafo. Es decir, si A y B están unidos por una arista y A está contaminado, la probabilidad de que B también lo esté es superior a la normal.

Se me ocurre probar esa hipótesis así:

library(igraph)

# mi grafo
g <- erdos.renyi.game(10000,
  p.or.m = 0.001, type="gnp")

min.mean.dist <- function(n){
  # contaminación al azar
  contaminados <- sample(V(g), n)

  # distancias entre aristas contaminadas
  res <- shortest.paths(g,
    v = contaminados, to = contaminados)
  diag(res) <- Inf

  # distancia al contaminado más próximo
  min.dist <- apply(res, 1, min, na.rm = T)

  # y su media
  mean(min.dist)
}

# histograma bajo la hipótesis nula
res <- replicate(100, min.mean.dist(100))

El resto son detalles que el lector atento sabrá completar por su cuenta.

Por si alguien aún no lo sabe: estamos todos citados en Albacete los días 17 y 18 de noviembre de 2016 en las VIII Jornadas de Usuarios de R.

Los detalles, aquí.

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.

En estadística, una muestra representativa tiene que contener las características relevantes de la población en las mismas proporciones en que están incluidas en tal población (referencia).

En estos tiempos, se están poniendo de moda alternativas a la muy tradicional democracia representativa que, en contraposición a ella, no aspiran a serlo. Y su principal problema radica, precisamente, en que no lo son.

Lo anterior no es más que una opinión: es la constatación de un hecho. Esta semana pasada, en aras de una versión más directa y asamblearia de la democracia, ha habido en mi barrio un par de eventos en los que en presencia de la alcaldesa de Madrid el uno y del concejal de mi distrito el otro, se han tratado temas que me interesan directamente. Pero, oh, fatalidad, a la hora en que yo (y muchos otros) estamos lejos y ocupados ganándonos el pan.

Muchas cosas han cambiado en los últimos 20 años. De hecho, acaba de hacer 20 años desde la primera vez que me conecté a internet y que tuve una cuenta de correo electrónico.

Sin embargo,

(Las cifras anteriores corresponden a taxis de las capitales de provincia; hay datos para otras en el INE).

Incidentalmente, hoy se han manifestado taxistas de toda España en Madrid para exigir protección para el sector.

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.

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