Mezclas

Estamos acostumbrados a la caracterización habitual de la distribución binomial negativa como el aburrido número de fracasos en una serie de ensayos de Bernoulli hasta lograr $r$ éxitos. Esto, junto con un poco de matemáticas de primero de BUP —todo aquello de combinaciones, etc.— lleva a la expresión conocida de su función de probabilidad,

$$\binom{n + x - 1}{x} p^r (1 - p)^x.$$

Pero esta caracterización, muy útil para resolver problemas de probabilidad construidos artificialmente para demostrar que los alumnos han estudiado la lección con aprovechamiento, se queda muy corta a la hora de proporcionar intuiciones sobre cómo, cuándo y por qué utilizarla en el ámbito en el que es más útil: el análisis de los procesos puntuales.

Cuando mezclas agua y tierra obtienes barro, una sustancia que comparte propiedades de sus ingredientes. Eso lo tenía muy claro de pequeño. Lo que en esa época me sorprendió mucho es que el agua fuese una mezcla de oxígeno e hidrógeno: ¡era muy distinta de sus componentes!

Porque no era una mezcla, obviamente. Era una combinación. En una combinación emergen propiedades inesperadas. Las mezclas, sin embargo, son más previsibles.

Pensaba en esto mientras escribía sobre la regularización de modelos (ridge, lasso y todas esas cosas). La regularización puede interpretarse como una mezcla de dos modelos: el original y el nulo (con todos los coeficientes iguales a cero). El modelo original tiene poco sesgo y mucha varianza; el nulo, prácticamente nada de varianza y muchísimo sesgo. El regularizado queda a medio camino. El original tiene varios, tal vez muchos, grados de libertad mientras que el nulo, ninguno (¿o uno?); puede considerarse que el número de grados de libertad del regularizado queda a medio camino.

A primeros de julio impartí un curso de estadística bayesiana aplicada con Stan. Tengo que examinar a los alumnos y he aquí el primero de los ejercicios:

En un país, se extrae una muestra de 2000 hombres y mujeres con la siguiente distribución:

men   <- 170 + 3 * rt(1000, 6)
women <- 160 + 2 * rt(1000, 5)
heights <- c(men, women)

Ajusta una distribución (una mezcla de dos distribuciones de Student) usando los datos anteriores, i.e., heights. Puedes suponer conocidos:

• Los pesos de la mezcla (0.5) cada uno.
• Que los grados de libertad de las t’s están entre 3 y 8 aproximadamente.
• Experimenta con otros tamaños muestrales y comenta los resultados obtenidos (y los tiempos de ejecución).

Nota: este problema está motivado por una aplicación real: el ajuste de distribuciones de pérdida en banca y seguros. Típicamente, se mezclan dos distribuciones, una para la cola de la distribución y otra para el cuerpo. Hay técnicas frecuentistas (p.e., EM) para resolver estos problemas. Pero me parecen menos naturales y menos flexibles que la ruta 100% bayesiana.

Pues no se sabe bien. Además, habrá quién pudiéndola haber averiguado, prefirió dejarse llevar por la intuición y errar. Pero volvamos a los hechos. Dado

En un país hipotético, las familias tienen críos hasta que nace el primer varón. En un año, en promedio, nacen:

— Carlos Gil Bellosta (@gilbellosta) December 10, 2017

la pregunta urgente es: ¿cuántos podrían haber conocido la respuesta? Suponiendo que el conocimiento de la respuesta es algo binarizable (¿lo es?), la distribución del número de respuestas correctas sería $latex pN + X$, donde $latex N$ es el número total de respuestas, $latex p$ es la proporción de quienes sabe la respuesta y $latex X \sim B(N - pN, 1/3)$, suponiendo siempre que $latex pN$ es entero.

[Tiempo después de la publicación de esta entrada hice otra, esta, en la que se ahonda en la función de pérdida usada en la reconstrucción del estilo o textura de las imágenes y que en esta serie no se trató con el detalle que el asunto requiere.]

En esta tercera entrada de la serie (aquí está la primera y la segunda) quiero ocuparme de las que llamé $latex f_1$ y $f_2$, las funciones involucradas. Que son las que obran la magia, por supuesto. Con casi cualquier otra opción se habría obtenido una patochada, pero estas son funciones especiales.

Siguiendo con el tema de la entrada de ayer, voy a tomar un vector $latex x_1$ tal como

y un vector $latex x_2$ como, por ejemplo,

para, con el concurso de unas funciones que revelaré mañana, obtener la siguiente mezcla de ambos:

Pas mal!

Arranco con esta una serie que estimo que será de tres entradas sobre cómo mezclar vectores con una aplicacioncilla que tal vez sorprenda a alguno.

Comenzaré fijando un vector $latex x_1 \in R^n$ y una función casi biyectiva $latex f_1:R^n \mapsto R^m$ todo lo suave (continua, diferenciable, etc.) que nos dé la gana. Casi no es un concepto matemático; el concepto propiamente matemático usaría el prefijo cuasi-, pero espero que se me permita seguir y prometo que lo que quiero dar a entender quedará claro más adelante.

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