Poisson

Fueron mis modelos favoritos un tiempo, cuando modelaba visitas y revisitas de usuarios a cierto malhadado portal.

Si las visitas fuesen aleatorias (en cierto sentido), tendrían un aspecto no muy distinto del que se obtiene haciendo

library(IHSEP)

suppressWarnings(set.seed(exp(pi * complex(imaginary = 1))))

tms <- simPois(int = function(x) .1, cens = 1000)
hist(tms, breaks = 100, main = "Proceso homogéneo de Poisson",
      xlab = "", ylab = "frecuencia")

Es decir,

o bien una distribución uniforme en el tiempo. Pero bien puede ocurrir que una visita incremente la probabilidad de otra inmediatamente después, por lo que las visitas tenderían a arracimarse en determinados momentos. Con el paquete [IHSEP](https://cran.r-project.org/package=IHSEP) de R pueden simularse (y ajustarse) este tipo de modelos. Por ejemplo,

Esta entrada muestra cómo afrontar (con Stan) un problema que encontré el otro día en un lugar que no puedo mencionar pero en el que sé que me leen (y los destinatarios sabrán que va por ellos).

El contexto es el siguiente: se hace un test A/B donde la variable de interés son unos conteos. Hay varios grupos (aquí los reduciré a dos) y los datos siguen aproximadamente (aquí omitiré la parte de la inflación de ceros) una distribución de Poisson. Pero solo aproximadamente: existe sobredispersión, es decir, la varianza de los datos excede su media.

Escribo esta entrada con cierta prevención porque soy consciente de que dan pábulo a determinadas teorías conspiranoicas de las que soy declarado enemigo. Pero es que los números de muertos en carretera por accidente en España en los últimos años,

(extraídos de aquí) dan que pensar: la varianza de las observaciones correspondientes a los años 2013, 2014 y 2015 es muy baja, demasiado baja. Al menos, si se da como bueno un modelo de Poisson para modelar esos conteos.

Tenemos un dato y un valor de referencia. Por ejemplo, el valor predicho por uno modelo y el observado. Queremos medir la distancia entre ambos. ¿En qué unidades?

Antes de eso, incluso, ¿para qué queremos medir esa distancia? Esta es la pregunta fácil: para ver cómo encaja en el modelo propuesto, para ver cómo lo sorprende, para cuantificar la perplejidad.

Los estadísticos están acostumbrados a medir la perplejidad en unas unidades que solo ellos entienden, si es que las entienden: desviaciones estándar. El z-score de un residuo es el número de desviaciones estándar que lo separan de su estimación. Si es una, exclaman ¡bah!; si es dos, ¡oh!; si es tres, ¡oooh!; si es cuatro, ¡ooooooh, válgame Dios!, etc.

Imagínate que quieres estabilizar la varianza (¡para qué!) de una distribución de Poisson. Los libros viejunos te dirán que saques la raíz cuadrada de tus valores.

Si en lugar de mirar en libros viejunos prestas atención a tus propios ojos, harás algo parecido a:

lambdas <- -10:10
lambdas <- 2^lambdas
res <- sapply(lambdas,
    function(lambda) sd(sqrt(rpois(1e5, lambda))))

para obtener

y averiguar dónde funciona y dónde no.

Si usas la transformación $latex f(x) = x^{2/3}$, como recomiendan en cierto artículo que no viene a cuento identificar, harás

En efecto,

mean(rpois(100000, 28 * 60 / 365) >= 10)
#[1] 0.01964

Por referencia,

• 28 es el número de días de febrero
• 60 viene de aquí
• 10 viene de aquí

Estimados señores:
Llevo 10 años revisando sus "CAJAS DE 100 CERILLAS"
En 3409 ocasiones he contado 99 o 101

😨 ¿ESTÁN USTEDES LOCOS? 😠 pic.twitter.com/hyqI9Ncxqg

— ☢️ 𝙍𝙖𝙙𝙞𝙖𝙘𝙩𝙞𝙫𝙤𝙈𝙖𝙣 ☢️ (@RadiactivoMan) February 16, 2017

Esta entrada, obviamente, viene a cuento de esta otra.

La distribución binomial (de parámetro n, p) es una suma de n variables aleatorias de Bernoulli independientes de parámetro p.

Independientes, reitero.

La distribución de Poisson es aproximadamente, una distribución binomial con un n muy grande y un p muy pequeño.

Los eventos subyacentes siguen siendo independientes, reitero.

Viene esto al caso de una tabla que ha circulado por Twitter,

en la que se comparan estimaciones de los parámetros $latex \lambda$ de una serie de distribuciones de Poisson… como si todas lo fuesen.

La distribución de Poisson se utiliza de oficio cuando se quiere modelar datos relativos a conteos. Sin embargo, tiene un problema serio: la varianza está fijada a la media: ambas son $latex \lambda$, el parámetro de la distribución.

Muy frecuentemente se observan datos con sobredispersión. Si $latex \lambda$ es 1000, el número esperado de eventos está contenido en un intervalo demasiado estrecho,

qpois(c(0.025, 0.975), 1000)
#[1]  938 1062

como para ser realista en muchas aplicaciones.

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