Estadística Bayesiana

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, $n_i$, que realiza el $i$-ésimo portero tiene una distribución binomial

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

donde $N_i$ es el número de disparos entre los palos y $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 $y_{1i}$ e $y_{2i}$ donde

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

e interesa saber si $\mu_1 = \mu_2$. Obviamente, se desconoce $\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 $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.

Hoy voy a hablar de esa especie de oxímoron que es el bootstrap bayesiano. Comenzaré planteando un pequeño problema bien conocido: tenemos números $x_1, \dots, x_n$ y hemos calculado su media. Pero nos preguntamos cómo podría variar dicha media (de realizarse otras muestras).

La respuesta de Efron (1979) es esta:

replicate(n, mean(sample(x, length(x), replace = TRUE)))

Es decir, crear muestras de $x_i$ con reemplazamiento y hacer la media de cada una de ellas para obtener su presunta distribución (o una muestra de la presunta distribución de esa media).

Te encargan un modelo. Por ejemplo, relacionado con el uso de tarjetas de débito y crédito (aunque a lo que me referiré ocurre en mil otros contextos). Una variable que consideras importante es la proporción de veces que se usa para sacar dinero de cajeros (y no para pagar en establecimientos). Así que, para cada cliente, divides el número de retiradas por el número de veces que la tarjeta se ha usado y obtienes ese número entre el 0 y el 1 (o entre el 0% y el 100%).

Una red bayesiana es algo de lo que ya hablé (y que me está volviendo a interesar mucho últimamente). En esencia, es un modelo probabilístico construido sobre un grafo dirigido acíclico.

Que, a su vez, es algo parecido a

que es un grafo (obviamente), dirigido (tiene flechas) y acíclico porque siguiéndolas no se llega nunca al punto de partida. Se puede construir modelos probabilísticos sobre ellos. Basta con definir para cada nodo $x$ la probabilidad condicional $P(x|A(x))$, donde $A(x)$ son sus padres directos. Con estas probabilidades condicionales (y un poco de esfuerzo) se puede construir la función de probabilidad completa, $P(x_1, \dots, x_n)$.