Estadística Bayesiana

Un banco tiene clientes. Los clientes usan la tarjeta de débito. La pueden usar de dos maneras: en cajero o para pagar (por productos y servicios). De cada cliente se tiene una secuencia de transacciones, etiquetadas como 1 o 0 según la use en cajero o no.

Para cada cliente, la secuencia de transacciones (más o menos larga) puede considerarse una secuencia intercambiable y, de acuerdo con el teorema de representación de de Finetti,

Mi búsqueda de ejemplos de aplicaciones con prioris informativas me ha conducido a Physiological pharmacokinetic analysis using population modeling and informative prior distributions, un artículo en el que se plantea un modelo jerárquico con dos tipos de distribuciones a priori:

Distribuciones muy informativas. Por ejemplo, el parámetro que representa la proporción del peso del hígado en un adulto, alrededor del 3.3% en promedio, que se modela con una distribución centrada en ese valor y una desviación estándar baja.

La estadística bayesiana se enseña en cursos de estadística (y, frecuentemente, envuelto en un aparataje matemático tan ofuscante como innecesario). Lo malo es que en los cursos y textos de estadística no existe información previa. La información previa sobre los fenómenos en los que se utilizaría la estadística bayesiana están en las aplicaciones, extramuros del muy agnóstico mundo de la estadística y la matemática.

Por eso, a los autores de los libros de estadística bayesiana y quienes enseñan cursos sobre lo mismo, enfrentados al problema de llenar de sentido la problemática distribución a priori, no se les ocurre nada mejor que discutir muy sesudamente la excepción (la priori no informativa) en lugar de la regla (la priori informativa). Reto al lector escéptico a que repase cualquier manual en la materia (que no haya sido escrito por Gelman) y compare el espacio que dedican a la selección de prioris no informativas con el de convenir una priori informativa decente.

La búsqueda de la causa más probable de un efecto tiene un nombre: razonamiento abductivo. Que el visitante al enlace anterior aprenderá distinto del deductivo y el inductivo.

Y que los viejos de estas páginas reconocerán en esta entrada que la formaliza y cuantifica en un caso concreto.

Dentro de unos días voy a hablar de estadística bayesiana en Machine Learning Spain. Plantearé una distribución a priori muy poco informativa:

alfa ~ gamma(10, 1);
beta ~ gamma(10, 1);

Me estoy preparando sicológicamente para que alguien me dé guerrita con lo de la subjetividad de las distribuciones a priori. Si tal es el caso, replicaré lo que sigue.

Hace unos días quise replicar el análisis. Pero la URL de la que bajo los datos dejó de contener los de la liga del año anterior y cargó los correspondientes al inicio (¿dos jornadas? ¿tres?) de la actual. ¡Apenas había datos!

Vuelvo a lo de Casillas inspirándome en el primer ejemplo de este artículo de Gelman et al.

El planteamiento es el siguiente: el número de paradas, $latex n_i$ que realiza el $latex i$-ésimo portero tiene una distribución binomial

$$ n_i \sim B(N_i, p_i)$$

donde $latex N_i$ es el número de disparos entre los palos y $latex p_i$ es la habilidad innata del portero. Estas habilidades innatas siguen una distribución dada, la de habilidades innatas de los porteros de primera división, que podemos suponer que sigue una distribución beta

El habitual problema de la diferencia de medias suele formularse de la siguiente manera: hay observaciones $latex y_{1i}$ e $latex y_{2i}$ donde

$$ y_{ji} \sim N(\mu_j, \sigma)$$

e interesa saber si $latex \mu_1 = \mu_2$. Obviamente, se desconoce $latex \sigma$. De cómo resolvió Gosset el problema están los libros de estadística llenos. En R,

set.seed(1234)
N1 <- 50
N2 <- 50
mu1 <- 1
mu2 <- -0.5
sig1 <- 1
sig2 <- 1
 
y1 <- rnorm(N1, mu1, sig1)
y2 <- rnorm(N2, mu2, sig2)
 
t.test(y1, y2)
# Welch Two Sample t-test
#
# data:  y1 and y2
# t = 4.7059, df = 95.633, p-value = 8.522e-06
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
#   0.5246427 1.2901923
# sample estimates:
#   mean of x  mean of y
# 0.5469470 -0.3604705

En rstan,

A partir de los comentarios de Olivier Núñez a mi entrada anterior casi homónima, se nos ha ocurrido a ambos de forma independiente y simultánea una manera alternativa de calcular el intervalo: minimizando su longitud.

a <- 3
b <- 5
alfa <- 0.05

# versión de la entrada anterior:
f <- function(x){
  (dbeta(x[2], a, b) - dbeta(x[1], a, b))^2 +
    (pbeta(x[2], a, b) - pbeta(x[1], a, b) -1 +  alfa)^2
}

res <- optim(c(a/(a+b), a/(a+b)), f)
res$par
#[1] 0.08052535 0.68463436

# nueva versión
f.alt <- function(x){
  qbeta(x+0.95, a, b) - qbeta(x, a, b)
}

res.alt <- optim(0.025, f.alt)
qbeta(c(res.alt$par, res.alt$par + 0.95), a, b)
#[1] 0.08054388 0.68464900

Tengo un parámetro, la p de una binomial, que supongo distribuido según una beta. Me da igual para el caso si la distribución a priori es o no informativa. Solo digo que la distribución a posteriori es otra beta con parámetros a y b.

Quiero construir un intervalo de credibilidad para p, es decir, encontrar un subintervalo de [0,1]

• dentro del cual la densidad de la beta sea mayor que fuera y que
• capture $latex 1-\alpha$ de la probabilidad total.

Gráficamente,

Esta entrada trata de cuadrados. Tales como estos

Son dos cuadrados de area 10 y 2.

En realidad, mi entrada trata de una configuración de cuadrados solo marginalmente más complicada, esta:

Todo el mundo podría decir (y es cierto) que el área de la intersección de los cuadrados es el 3.3% de la del mayor y el 16.5% de la del menor. Son dos afirmaciones ambas ciertas y, por supuesto, compatibles.