1 Modelo general
Si Su estado de información acerca de las secuencia de variables binarias \(y_1,\ldots,y_n\) es intercambiable, entonces el modelamiento \(y_1,\ldots,y_n\) admite representación jerárquica de la forma: \[\begin{align} y_i\mid\theta &\stackrel{\text{iid}}{\sim}\textsf{Ber}(\theta)\,,\quad i = 1,\ldots,n \\ \theta &\sim p(\theta) \end{align}\] Así, la distribución muestral (distribución condicional conjunta) de \(\boldsymbol{y} = (y_1,\ldots,y_n)\) dado \(\theta\) está dada por \[ p(\boldsymbol{y}\mid\theta) = \theta^{s}(1-\theta)^{n - s}\,, \] donde \(s = \sum_{i=1}^n y_i\), lo cual sugiere que \(s\) es un estadístico suficiente para \(\theta\) (i.e., \(s\) contiene toda la información que proviene de los datos acerca de \(\theta\); saber \(s\) es suficiente para hacer inferencia sobre \(\theta\)). La demostración formal se puede hacer por medio del Teorema de Factorización de Fisher-Neyman.
Por lo tanto, la distribución posterior es \[ p(\theta\mid\boldsymbol{y}) \propto \theta^{s}(1-\theta)^{n - s}\,p(\theta)\,, \]
Dado que las \(y_i\)’s son condicionalmente i.i.d. dado \(\theta\) y \(s\) es un estadístico suficiente para \(\theta\), entonces el modelo es equivalente a \[\begin{align*} s\mid\theta &\sim \textsf{Bin}(n,\theta) \\ \theta &\sim p(\theta) \end{align*}\]
2 Familias conjugadas
La idea básica consiste en encontrar una familia de distribuciones \(\mathcal{P}\) de tal forma que el producto de los miembros de esta familia con la distribución muestral también sea parte de \(\mathcal{P}\).
(Definición.) Una familia de distribuciones \(\mathcal{P}\) se denomina conjugada para una distribución muestral dada \(p(\boldsymbol{y}\mid\boldsymbol{\theta})\), si \(p(\boldsymbol{\theta}) \in \mathcal{P}\) entonces \(p(\boldsymbol{\theta}\mid \boldsymbol{y}) \in \mathcal{P}\).
Las previas conjugadas conllevan a cálculos fáciles de realizar, pero algunos casos pueden ser poco flexibles para representar Su estado de información previa.
La distribución previa no tiene que ser necesariamente conjugada.
(Propiedad.) Utilizando familias conjugadas, la media posterior se puede expresar como un promedio ponderado de la media previa y la media muestral con pesos proporcionales al tamaño de la muestra previa y el tamaño de la muestra.
(Teorema.) Dada una distribución muestral miembro de la familia exponencial, cualquier información previa se puede expresar como una mezcla de previas conjugadas (Diaconis y Ylvisaker 1985).
3 Modelo Beta-Binomial
La familia de distribuciones Beta es conjugada para la distribución muestral Binomial.
La variable aleatoria \(X\) tiene distribución Beta con parámetros \(\alpha,\beta > 0\), i.e., \(X\mid\alpha,\beta\sim\textsf{Beta}(\alpha,\beta)\), si su función de densidad de probabilidad es \[ p(x\mid\alpha,\beta) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\,\Gamma(\beta)}\,x^{\alpha-1}\,(1-x)^{\beta-1},\quad 0<x<1\,, \] donde \(\Gamma(x)=\int_0^\infty u^{x-1}\,e^{-u}\,\text{d}u\) es la función Gamma.
Así, el modelo Beta-Binomial es \[\begin{align*} s\mid\theta &\sim \textsf{Bin}(s\mid n,\theta) \\ \theta &\sim \textsf{Beta}(a,b) \end{align*}\] donde \(a\) y \(b\) son los hiperparámetros (cantidades fijas conocidas) del modelo. Los hiperparámetros se eligen de tal forma que la distribución previa refleje Su estado de información previo.
Bajo el modelo Beta-Binomial se tiene que la distribución posterior es \[ \theta\mid s \sim \textsf{Beta}(\theta\mid a + s, b+n-s)\,. \]
Además, la distribución marginal de \(s\) es \[ p(s) = \frac{\Gamma(n+1)}{\Gamma(s+1)\Gamma(n-s+1)}\,\frac{\Gamma(a+b)}{\Gamma(a)\,\Gamma(b)}\,\frac{\Gamma(a+s)\,\Gamma(b+n-s)}{\Gamma(a+b+n)}\,. \] Esta distribución se denomina distribución Beta-Binomial (mezcla de distribuciones Binomial, usando como pesos valores de una distribución Beta).
De otra parte, la media posterior es \[ \textsf{E}(\theta\mid s) = \frac{a+s}{a+b+n} = \frac{a+b}{a+b+n}\cdot\frac{a}{a+b}+\frac{n}{a+b+n}\cdot\frac{s}{n}\,, \] la cual es un promedio ponderado del valor esperado previo y la media muestral.
- Tal observación conlleva a la siguiente interpretación de los hiperparámetros:
- \(a\) = número previo de 1’s.
- \(b\) = número previo de 0’s.
- \(a+b\) = tamaño muestral previo.
- Si \(n >> a+b\), entonces la mayoría de la información sobre \(\theta\) proviene de los datos en lugar de la información previa.
(Propiedad.) Los resultados Bayesianos son aproximadamente iguales a los frecuentistas si: 1. la distribución previa es difusa (no informativa) y 2. el tamaño de la muestra es “grande”. En este caso está bien hacer un cálculo frecuentista e interpretarlo como uno Bayesiano.
Con tamaños de muestra pequeños, las propiedades frecuentistas no son buenas y la alternativa Bayesiana constituye una mejor alternativa.
Mal Bayesiano \(\leq\) Frecuentista \(\leq\) Buen Bayesiano (la previa es fuerte y en concordancia con el mundo).
También, la distribución predictiva posterior es tal que \[ \textsf{Pr}(y^*=1\mid s) = \frac{a+s}{a+b+n}\,. \] La distribución predictiva no depende de cantidades desconocidas. La distribución predictiva depende de los datos. Por tal motivo, \(y^*\) no es independiente de \(s\).
4 Regiones de confianza
(Definición.) Una región de confianza para \(\boldsymbol{\theta}\) es una región del espacio de parámetros \(\Theta\) que contiene a \(\boldsymbol{\theta}\) con alta probabilidad.
Por ejemplo, en el caso univariado, luego de observar los datos se construye un intervalo \((l,u)\) tal que \(\textsf{Pr}(l<\theta<u\mid \boldsymbol{y})\) sea alta.
Nivel de confianza Bayesiano: probabilidad de que el intervalo cubra el verdadero valor del parámetro después de que los datos sean observados.
Nivel de confianza Frecuentista: probabilidad de que el intervalo cubra el verdadero valor del parámetro antes de que los datos sean observados. ¿A qué es igual esta probabilidad una vez se han observado los datos?
En el caso univariado, la manera más sencilla de obtener intervalos credibilidad en el paradigma Bayesiano es por medio de los percentiles de la distribución posterior: \((\theta_{\alpha/2},\theta_{1-\alpha/2}\mid \boldsymbol{y})\) es una región de confianza al \(100(1-\alpha)\%\) basado en cuantiles, y por lo tanto \[ \textsf{Pr}\left(\theta_{\alpha/2} < \theta < \theta_{1-\alpha/2}\mid \boldsymbol{y}\right) = 1-\alpha\,. \]
(Teorema.) Un intervalo de credibilidad con un nivel de confianza Bayesiano de \((1-\alpha)\%\), también tiene asintóticamente un nivel de confianza Frecuentista de \((1-\alpha)\%\) (Hartigan 1966).
Es posible investigar las propiedades Frecuentistas de un intervalo de credibilidad por medio de simulación.
5 Ejemplo
A cada mujer de 65 años o más en la Encuesta Social General de 1998 en Estados Unidos se le preguntó si en general se consideraban felices o no. Se define \(y_i = 1\) si el individuo \(i\) informó sentirse generalmente feliz, y \(y_i = 0\) en caso contrario, \(i = 1,\ldots,n\).
# cargar data
Y <- read.csv("GSS1998.txt")
# extraer variable
y <- Y[(Y$YEAR == 1998) & (Y$AGE >= 65) & (Y$FEMALE == 1),]$HAPUNHAP
y[y > 4] <- NA
y[y <= 2] <- 1
y[y > 2] <- 0
y <- y[!is.na(y)]
# frecuencias
table(y)## y
## 0 1
## 11 118# tamaño de muestra
n <- length(y)
print(n)## [1] 129# estadistico suficiente
sy <- sum(y)
print(sy)## [1] 118# previa Beta(1,1)
a <- 1
b <- 1
# parametros de la posterior
ap <- a + sy
bp <- b + n - sy
# grafico
theta <- seq(0, 1, length = 1000)
plot(theta, dbeta(theta, ap, bp), type = "l", col = 4, xlab = expression(theta),
ylab = expression(paste("p","(",theta," | ",y,")",sep="")), xlim = c(0.8, 1))
lines(theta, rep(1, length(theta)), type = "l", col = 1)
legend("topright", legend = c("Previa", "Posterior"), col = c(1, 4), lty = 1,
lwd = 2, bty = "n")
# moda posterior
(ap-1)/(ap + bp - 2)## [1] 0.9147287# media posterior
ap/(ap + bp)## [1] 0.9083969# mediana posterior
(ap - 1/3)/(ap + bp - 2/3)## [1] 0.9104859# varianza posterior
(ap*bp)/((ap + bp)^2*(ap + bp + 1))## [1] 0.0006303934# intervalo de credibilidad al 95%
qbeta(c(.025,.975), shape1 = a+sy, shape2 = b+n-sy)## [1] 0.8536434 0.9513891# grafico
theta <- seq(0, 1, length = 1000)
plot(theta, dbeta(theta, ap, bp), type = "l", col = 4, xlab = expression(theta),
ylab = expression(paste("p","(",theta," | ",y,")",sep="")), xlim = c(0.8, 1))
lines(theta, rep(1, length(theta)), type = "l", col = 1)
abline(v = qbeta(c(.025,.975), a+sy,b+n-sy), lty = 2, lwd = 2, col = 2)
abline(v = ap/(ap + bp), lty = 2, lwd = 2, col = 3)
legend("topright", legend = c("Previa", "Posterior", "IC 95%", "Media"),
col = c(1, 4, 2, 3), lty = 1, lwd = 2, bty = "n")
# probabilidad a priori de que theta > 0.9
pbeta(q = 0.9, shape1 = a, shape2 = b, lower.tail = F)## [1] 0.1# probabilidad a posteriori de que theta > 0.9
pbeta(q = 0.9, shape1 = a+sy, shape2 = b+n-sy, lower.tail = F)## [1] 0.6576287