La regla de Bayes es el núcleo de la estadística bayesiana. Permite actualizar la probabilidad de una hipótesis a partir de la evidencia observada. En este capítulo se presenta el teorema de Bayes sobre sistemas completos de sucesos, se introducen los conceptos de distribución a priori, función de verosimilitud y distribución a posteriori, y se simula el proceso de actualización con R. Al final se estudian las distribuciones conjugadas, que simplifican mucho el cálculo de la posterior.
Definición 4.1 (Sistema completo de sucesos) Una colección de sucesos \(A_1,A_2,\ldots,A_n\) de un mismo espacio muestral \(\Omega\) es un sistema completo si cumple las siguientes condiciones:
En realidad, un sistema completo de sucesos es una partición del espacio muestral de acuerdo con algún atributo, como el sexo o el grupo sanguíneo.
Conocer las probabilidades de un determinado suceso en cada una de las partes de un sistema completo puede ser útil para calcular su probabilidad global.
Teorema 4.1 (Probabilidad total) Dado un sistema completo de sucesos \(A_1,\ldots,A_n\) de un espacio muestral \(\Omega\) y un suceso \(B\) de ese mismo espacio muestral, la probabilidad de \(B\) se puede calcular mediante la fórmula
\[P(B) = \sum_{i=1}^n P(A_i\cap B) = \sum_{i=1}^n P(A_i)P(B|A_i).\]
Ejemplo 4.1 Un test diagnóstico para una enfermedad \(E\) tiene una sensibilidad del 90 % y una especificidad del 95 %, es decir,
\[ P(+|E)=0{,}9 \quad \mbox{y} \quad P(-|\overline{E})=0{,}95. \]
Si la prevalencia de la enfermedad en la población es del 10 %, ¿qué porcentaje de positivos se espera en una muestra aleatoria de la población?
La siguiente tabla de probabilidades conjuntas resume la situación:
| Enfermo | No enfermo | Total | |
|---|---|---|---|
| \(+\) | 0,090 | 0,045 | 0,135 |
| \(-\) | 0,010 | 0,855 | 0,865 |
| Total | 0,100 | 0,900 | 1 |
Para responder a la pregunta se aplica el teorema de la probabilidad total con el sistema completo \(\{E,\overline{E}\}\):
\[ P(+) = P(E)P(+|E)+P(\overline E)P(+|\overline E) = 0{,}1\cdot 0{,}9 + 0{,}9\cdot 0{,}05 = 0{,}135. \]
Es decir, alrededor del 13,5 % de los resultados del test serán positivos.
En el fondo se trata de una media ponderada de probabilidades.
Los sucesos de un sistema completo \(A_1,\ldots,A_n\) también pueden verse como las distintas hipótesis ante un determinado hecho \(B\). En estas condiciones resulta útil poder calcular las probabilidades a posteriori \(P(A_i|B)\) de cada una de las hipótesis.
Teorema 4.2 (Teorema de Bayes) Dado un sistema completo de sucesos \(A_1,\ldots,A_n\) y un suceso \(B\) de un mismo espacio muestral \(\Omega\), la probabilidad de cada suceso \(A_i\) condicionada por \(B\) puede calcularse con la fórmula
\[P(A_i|B) = \frac{P(A_i\cap B)}{P(B)} = \frac{P(A_i)P(B|A_i)}{\sum_{j=1}^n P(A_j)P(B|A_j)}.\]
Ejemplo 4.2 En el ejemplo anterior, una pregunta más interesante es qué diagnosticar a una persona que ha dado positivo en el test. Podemos interpretar \(E\) y \(\overline{E}\) como las dos posibles hipótesis ante el resultado positivo del test.
Probabilidades a priori (sin información sobre el resultado del test)
\[ P(E)=0{,}1 \quad \mbox{y}\quad P(\overline E)=0{,}9. \]
Sin información, el diagnóstico sería que la persona no tiene la enfermedad.
Sin embargo, si al realizar el test a una persona se obtiene un resultado positivo, esa información condiciona las hipótesis. Para decidir entre ellas hay que calcular sus probabilidades a posteriori, \(P(E|+)\) y \(P(\overline{E}|+)\):
\[ \begin{align*} P(E|+) &= \frac{P(E)P(+|E)}{P(E)P(+|E)+P(\overline{E})P(+|\overline{E})} = \frac{0{,}1\cdot 0{,}9}{0{,}1\cdot 0{,}9 + 0{,}9\cdot 0{,}05} = \frac{0{,}09}{0{,}135}=0{,}667,\\ P(\overline{E}|+) &= \frac{P(\overline{E})P(+|\overline{E})}{P(E)P(+|E)+P(\overline{E})P(+|\overline{E})} = \frac{0{,}9\cdot 0{,}05}{0{,}1\cdot 0{,}9 + 0{,}9\cdot 0{,}05} = \frac{0{,}045}{0{,}135}=0{,}333. \end{align*} \]
Con un resultado positivo, la probabilidad de tener la enfermedad ha aumentado de 0,10 a 0,667, y el diagnóstico pasa a ser que la persona tiene la enfermedad.
Ejercicio 4.1 Se dispone de dos test, \(A\) y \(B\), para diagnosticar una enfermedad con una prevalencia del 1 % en la población. El test \(A\) tiene una sensibilidad del 98 % y una especificidad del 80 %, mientras que el test \(B\) tiene una sensibilidad del 75 % y una especificidad del 99 %.
Para cribado se prioriza una alta sensibilidad, es decir, detectar a todos los enfermos minimizando los falsos negativos. Por tanto, el test \(A\) es mejor para cribado. Su baja especificidad (80 %) produce falsos positivos, pero estos se filtran con una prueba confirmatoria posterior.
Aplicando el teorema de Bayes al resultado positivo del test \(A\):
\[ \begin{align*} P(E \mid +_A) &= \frac{P(+_A \mid E)\cdot P(E)}{P(+_A \mid E)\cdot P(E) + P(+_A \mid \overline{E})\cdot P(\overline{E})} \\ &= \frac{0{,}98 \times 0{,}01}{0{,}98 \times 0{,}01 + 0{,}20 \times 0{,}99} = \frac{0{,}0098}{0{,}2078} \approx 0{,}047 \end{align*} \]
Como \(P(E|+_A) = 0{,}047 < 0{,}5\), diagnosticaríamos sano, aunque derivaríamos a la persona a una prueba confirmatoria.
Ahora usamos \(P(E|+_A)=0{,}047\) en lugar de \(P(E)=0{,}01\) como nueva probabilidad a priori, y aplicamos Bayes con el test \(B\):
\[ \begin{align*} P(E \mid +_A, +_B) &= \frac{P(+_B \mid E)\cdot P(E \mid +_A)}{P(+_B \mid E)\cdot P(E \mid +_A) + P(+_B \mid \overline{E})\cdot P(\overline{E} \mid +_A)} \\ &= \frac{0{,}75 \times 0{,}047}{0{,}75 \times 0{,}047 + 0{,}01 \times 0{,}953} \approx 0{,}788. \end{align*} \]
Como \(P(E|+_A,+_B) \approx 0{,}79 > 0{,}5\), diagnosticaríamos enfermedad.
En terminología bayesiana, si \(H\) representa una hipótesis y \(D\) los datos observados, el teorema de Bayes se puede expresar como una regla de actualización de la probabilidad de la hipótesis \(H\) a la luz de la evidencia \(D\):
\[ \underbrace{P(H|D)}_{\text{Probabilidad a posteriori}} = \frac{\overbrace{P(H)}^{\text{Probabilidad a priori}}\times \overbrace{P(D|H)}^{\text{Verosimilitud}}}{\underbrace{P(D)}_{\text{Evidencia marginal}}} \]
El denominador \(P(D)\) actúa como normalizador: garantiza que las probabilidades a posteriori de todas las hipótesis sumen 1.
En la regla de Bayes, la probabilidad \(P(D|H)\) se conoce como verosimilitud de la hipótesis \(H\) dada la evidencia \(D\).
Cuando tenemos varias hipótesis, la verosimilitud de cada una de ellas se puede comparar con la de las demás para determinar cuál es la más compatible con los datos observados.
Definición 4.2 (Función de verosimilitud) Dado un sistema completo de sucesos con hipótesis \(H_1, H_2, \ldots, H_n\) y un suceso \(D\) que representa la evidencia observada, la función de verosimilitud asigna a cada hipótesis \(H_i\) la probabilidad de observar la evidencia \(D\) si esa hipótesis fuera cierta:
\[ L(H_i|D) = P(D|H_i). \]
Los valores de la función de verosimilitud son probabilidades condicionadas, pero la variable de la función de verosimilitud es el condicionante, es decir, la hipótesis \(H_i\), y no la evidencia \(D\):
\[ L(x|D) = P(D|x). \]
Ejemplo 4.3 Consideremos de nuevo el test diagnóstico con sensibilidad del 90 % y especificidad del 95 %. Dada la evidencia de que el resultado del test ha sido positivo, la función de verosimilitud de las dos hipótesis \(E\) y \(\overline{E}\) es
\[ L(E|+) = P(+|E) = 0{,}9 \quad \mbox{y} \quad L(\overline{E}|+) = P(+|\overline{E}) = 1 - 0{,}95 = 0{,}05. \]
En este caso, la hipótesis \(E\) es mucho más compatible con la evidencia \(+\) que la hipótesis \(\overline{E}\).
La función de verosimilitud se puede interpretar como una medida de la compatibilidad de cada hipótesis con los datos observados: cuanto mayor sea la verosimilitud de una hipótesis dada la evidencia, más plausible será la hipótesis.
Es importante destacar que la función de verosimilitud no es una función de probabilidad: la suma de las verosimilitudes de todas las hipótesis no tiene por qué ser igual a 1.
Ejemplo 4.4 La siguiente tabla recoge la probabilidad a priori y la función de verosimilitud de las hipótesis \(E\) y \(\overline{E}\) dada la evidencia \(+\) en el ejemplo del test diagnóstico.
| Hipótesis | Probabilidad a priori \(P(H)\) | Función de verosimilitud \(L(H\mid +)\) |
|---|---|---|
| \(E\) | 0,1 | 0,9 |
| \(\overline{E}\) | 0,9 | 0,05 |
Si disponemos de la probabilidad a priori de cada hipótesis y de su función de verosimilitud dada la evidencia, podemos calcular la probabilidad a posteriori de cada hipótesis con la regla de Bayes:
\[ P(H_i|D) = \frac{P(H_i)\times P(D|H_i)}{P(D)}. \]
La probabilidad del denominador \(P(D)\) es un factor común para todas las hipótesis, por lo que la probabilidad a posteriori de cada hipótesis es proporcional al producto de su probabilidad a priori por su verosimilitud:
\[ P(H_i|D) \propto P(H_i)\times P(D|H_i). \]
Ejemplo 4.5 En el ejemplo del test diagnóstico, dado que el test ha resultado positivo, obtenemos las siguientes plausibilidades a posteriori:
| Hipótesis | Probabilidad a priori \(P(H)\) | Verosimilitud \(L(H\mid +)\) | Plausibilidad a posteriori \(P(H)\times L(H\mid +)\) |
|---|---|---|---|
| \(E\) | 0,1 | 0,9 | 0,090 |
| \(\overline{E}\) | 0,9 | 0,05 | 0,045 |
La suma de las plausibilidades a posteriori no es igual a 1, por lo que no se pueden interpretar como probabilidades. Sin embargo, sí podemos compararlas: la hipótesis \(E\) es el doble de plausible que la hipótesis \(\overline{E}\) dada la evidencia de que el test ha resultado positivo.
¿Cómo podemos transformar las plausibilidades a posteriori en probabilidades a posteriori? Normalizándolas: dividiendo cada una de ellas por la suma de todas las plausibilidades.
| Hipótesis | Plausibilidad a posteriori | Probabilidad a posteriori \(P(H\mid +)\) |
|---|---|---|
| \(E\) | 0,090 | 0,090/(0,090+0,045) = 0,667 |
| \(\overline{E}\) | 0,045 | 0,045/(0,090+0,045) = 0,333 |
Podemos simular el procedimiento de actualización de la probabilidad de una hipótesis a la luz de la evidencia que realiza la regla de Bayes. Para ello necesitamos tomar muestras aleatorias de:
Al filtrar después la muestra conjunta por la evidencia observada, obtenemos una muestra de la distribución a posteriori.
Ejercicio 4.2 Obtener una muestra aleatoria de tamaño \(n=10000\) de la distribución a priori de la hipótesis \(H\), con probabilidades \(P(E)=0{,}1\) y \(P(\overline{E})=0{,}9\), y representar gráficamente la distribución de frecuencias muestral.
df_H = tibble(hipótesis = c("E", "NoE")) # Hipótesis
prior <- c(0.1, 0.9) # P(E) y P(No E)
set.seed(123) # Para reproducibilidad
n <- 10000 # Tamaño de la muestra
# Muestreo de la distribución a priori
muestra_prior <- slice_sample(df_H, n = n, weight_by = prior, replace = TRUE)
# Diagrama de barras de la distribución de frecuencias
ggplot(muestra_prior, aes(x = hipótesis)) +
stat_count(fill = myblue) +
labs(x = "Hipótesis", y = "Frecuencia absoluta") +
theme_minimal()
Para simular la distribución conjunta de la hipótesis y el resultado del test, generamos para cada individuo de la muestra un resultado del test según la función de verosimilitud: \(P(+|E)=0{,}9\) y \(P(+|\overline E)=0{,}05\).
# Muestreo de la función de verosimilitud (resultado del test)
muestra_prior <- muestra_prior |>
mutate(test = ifelse(runif(n()) < ifelse(hipótesis == "E", 0.9, 0.05), "+", "-"))
# Diagrama de barras de la distribución de frecuencias del test
ggplot(muestra_prior, aes(x = hipótesis, fill = test)) +
stat_count(position = "stack") +
labs(x = "Hipótesis", y = "Frecuencia absoluta") +
theme_minimal()
Una vez tenemos la distribución conjunta, podemos obtener una muestra de la distribución a posteriori filtrando los casos en los que el resultado del test ha sido positivo.
# Filtrar casos con test positivo
muestra_posterior <- muestra_prior |> filter(test == "+")
# Diagrama de barras de la distribución de frecuencias de la muestra a posteriori
ggplot(muestra_posterior, aes(x = hipótesis, y = after_stat(prop), group = 1)) +
stat_count(fill = myblue) +
labs(x = "Hipótesis", y = "Frecuencia relativa") +
theme_minimal()
Ejercicio 4.3 Repetimos la simulación con 100000 individuos. La prevalencia de la enfermedad es del 1 % (prior \(P(E)=0{,}01\)), el test \(A\) tiene sensibilidad \(0{,}98\) y especificidad \(0{,}80\), y el test \(B\) tiene sensibilidad \(0{,}75\) y especificidad \(0{,}99\).
df_H = tibble(hipótesis = c("Enfermo", "Sano")) # Hipótesis
prior <- c(0.01, 0.99) # P(Enfermo) y P(Sano)
set.seed(123) # Para reproducibilidad
n <- 100000 # Tamaño de la muestra
# Muestreo de la distribución a priori
muestra_prior <- slice_sample(df_H, n = n, weight_by = prior, replace = TRUE)
muestra_prior |> ggplot(aes(x = hipótesis)) +
stat_count(fill = myblue) +
labs(title = "Distribución a priori", x = "Hipótesis", y = "Frecuencia absoluta")
# Muestreo de la distribución conjunta con el test A
muestra_prior <- muestra_prior |>
mutate(test = ifelse(runif(n()) < ifelse(hipótesis == "Enfermo", 0.98, 0.2), "+", "-"))
ggplot(muestra_prior, aes(x = hipótesis, fill = test)) +
stat_count(position = "stack") +
labs(title = "Distribución de frecuencias del resultado del test A",
x = "Hipótesis", y = "Frecuencia absoluta") +
theme_minimal()
# Filtrar casos con test A positivo
muestra_posterior <- muestra_prior |> filter(test == "+")
ggplot(muestra_posterior, aes(x = hipótesis, y = after_stat(prop), group = 1)) +
stat_count(fill = myblue) +
labs(title = "Distribución a posteriori tras un test A positivo",
x = "Hipótesis", y = "Frecuencia relativa") +
theme_minimal()
# Muestreo de la distribución conjunta con el test B sobre los positivos de A
muestra_prior2 <- muestra_posterior |>
mutate(test = ifelse(runif(n()) < ifelse(hipótesis == "Enfermo", 0.75, 0.01), "+", "-"))
ggplot(muestra_prior2, aes(x = hipótesis, fill = test)) +
stat_count(position = "stack") +
labs(title = "Distribución de frecuencias del resultado del test B",
x = "Hipótesis", y = "Frecuencia absoluta") +
theme_minimal()
# Filtrar casos con test B positivo
muestra_posterior <- muestra_prior2 |> filter(test == "+")
ggplot(muestra_posterior, aes(x = hipótesis, y = after_stat(prop), group = 1)) +
stat_count(fill = myblue) +
labs(title = "Distribución a posteriori tras dos test positivos",
x = "Hipótesis", y = "Frecuencia relativa") +
theme_minimal()
Hemos visto que, para obtener la distribución de probabilidad a posteriori de una hipótesis dada la evidencia, hay que combinar la distribución a priori de la hipótesis con la función de verosimilitud dada la evidencia, mediante la regla de Bayes.
Existen familias de distribuciones de probabilidad que, al combinarse con determinadas verosimilitudes, mantienen su forma funcional. Estas familias se denominan distribuciones conjugadas y son de gran utilidad en la inferencia bayesiana, porque la posterior se obtiene en forma cerrada sin necesidad de calcular la integral del denominador. Las parejas más comunes son:
Si \(p\sim \mbox{Beta}(\alpha, \beta)\) y \(X\sim B(n, p)\), entonces la distribución a posteriori de \(p\) dado \(X=k\) es también una distribución beta:
\[ p|X=k \sim \mbox{Beta}(\alpha+k, \beta+n-k). \]
Ejemplo 4.6 Supongamos que queremos estimar la proporción de fumadores en nuestra ciudad y que tenemos una distribución a priori de \(p\) dada por una distribución beta con parámetros \(\alpha=2\) y \(\beta=5\), es decir, \(p\sim \mbox{Beta}(2, 5)\). Si tomamos una muestra aleatoria de \(n=30\) personas y observamos que \(k=10\) de ellas son fumadoras, la distribución a posteriori de \(p\) dado \(X=10\) es
\[ p|X=10 \sim \mbox{Beta}(2+10, 5+30-10) = \mbox{Beta}(12, 25). \]
Si \(\lambda \sim \mbox{Gamma}(\alpha, \beta)\) y \(X\sim \mbox{P}(\lambda)\), entonces la distribución a posteriori de \(\lambda\) dado \(X=k\) es también una distribución gamma:
\[ \lambda|X=k \sim \mbox{Gamma}(\alpha+k, \beta+1). \]
Si se toman \(n\) muestras independientes de la misma distribución \(P(\lambda)\), la distribución a posteriori de \(\lambda\) dado \(X_1=x_1,\ldots,X_n=x_n\) es también una distribución gamma:
\[ \lambda|X_1=x_1,\ldots,X_n=x_n \sim \mbox{Gamma}\left(\alpha+\sum_{i=1}^n x_i, \beta+n\right). \]
Ejemplo 4.7 Supongamos que queremos estimar la tasa de nacimientos diarios en una ciudad. Tenemos una distribución a priori de \(\lambda\) dada por una distribución gamma con parámetros \(\alpha=2\) y \(\beta=5\), es decir, \(\lambda \sim \mbox{Gamma}(2, 5)\). Si observamos que en una muestra de 5 días se han producido \(X_1=10, X_2=8, X_3=9, X_4=11, X_5=7\) nacimientos, la distribución a posteriori de \(\lambda\) es
\[ \begin{align*} \lambda |X_1=10, X_2=8, X_3=9, X_4=11, X_5=7 &\sim \mbox{Gamma}(2+10+8+9+11+7, 5+5) \\ &= \mbox{Gamma}(47, 10). \end{align*} \]
Si \(\mu \sim N(\theta, \tau^2)\) y \(X\sim N(\mu, \sigma^2)\), entonces la distribución a posteriori de \(\mu\) dado \(X=x\) es también una distribución normal:
\[ \mu|X=x \sim N\left(\frac{\theta\sigma^2 + x\tau^2}{\sigma^2 + \tau^2}, \frac{\sigma^2\tau^2}{\sigma^2 + \tau^2}\right). \]
Si se toman \(n\) muestras independientes de la misma distribución \(N(\mu, \sigma^2)\), la distribución a posteriori de \(\mu\) dado \(X_1=x_1,\ldots,X_n=x_n\) es también normal:
\[ \mu|\{X_1=x_1,\ldots,X_n=x_n\} \sim N\left(\frac{\theta\sigma^2 + \bar{x}\,n\tau^2}{\sigma^2 + n\tau^2}, \frac{\sigma^2\tau^2}{\sigma^2 + n\tau^2}\right), \]
donde \(\bar{x}\) es la media muestral.
Ejemplo 4.8 Supongamos que el peso perdido con una dieta sigue una distribución normal \(N(\mu, \sigma^2=4)\) y queremos estimar el peso medio \(\mu\) perdido con la dieta. Tenemos una distribución a priori para \(\mu\) dada por una normal con media \(\theta=10\) y varianza \(\tau^2=1\), es decir, \(\mu \sim N(10, 1)\). Si tomamos una muestra aleatoria de \(n=20\) personas que han seguido la dieta y observamos que el peso medio perdido es \(\bar{x}=8\), la distribución a posteriori de \(\mu\) es
\[ \mu|\{X_1,\ldots,X_{20}\} \sim N\left(\frac{10\cdot 4 + 8\cdot 20\cdot 1}{4 + 20\cdot 1}, \frac{4\cdot 1}{4 + 20\cdot 1}\right) = N(8{,}33;\ 0{,}17). \]
Un modelo bayesiano es un modelo estadístico que incorpora la incertidumbre sobre los parámetros desconocidos mediante distribuciones de probabilidad a priori, y que utiliza la regla de Bayes para actualizar estas distribuciones a posteriori a medida que se obtiene nueva información a través de los datos observados.
Los modelos bayesianos pueden considerarse modelos de aprendizaje automático, porque aprenden distribuciones con la experiencia que se va adquiriendo a partir de los datos observados.
Para definir un modelo bayesiano es necesario especificar:
Dependiendo de la naturaleza de las hipótesis o parámetros desconocidos, el espacio de hipótesis puede ser discreto o continuo, lo que da lugar a dos tipos de distribuciones a priori:
Del mismo modo, según la naturaleza de los datos observados, la función de verosimilitud puede modelarse mediante distribuciones de probabilidad discretas o continuas:
Según el tipo de la distribución a priori y de la distribución de los datos que se usa para calcular la función de verosimilitud, tenemos distintos tipos de modelos bayesianos.
Cuando el parámetro de interés es una proporción \(\theta\), como la proporción de fumadores en una población, y los datos observados son el número de éxitos \(X\) en una serie de \(n\) ensayos independientes, como el número de personas fumadoras en una muestra, el modelo bayesiano más utilizado es el modelo beta-binomial:
Ejercicio 4.4 Supongamos que queremos estimar la proporción de fumadores en nuestra ciudad. En varias muestras tomadas en ciudades similares, la proporción media estaba alrededor del 33 %, y la menor y la mayor proporción observadas fueron del 19 % y del 48 %, respectivamente.
¿Qué distribución beta podríamos utilizar para modelizar la distribución a priori?
Puedes usar la aplicación distribuciones para probar distintas distribuciones beta y encontrar la que mejor se ajuste a la información.
Una distribución \(\mbox{Beta}(30, 60)\) tiene media \(30/90 \approx 0{,}33\) y se concentra en el rango de proporciones observado.
library(bayesrules)
plot_beta(alpha = 30, beta = 60) +
labs(x = "Proporción de fumadores (θ)", y = "Densidad")
Ejercicio 4.5 Tras realizar un nuevo sondeo en nuestra ciudad, hemos observado que en una muestra de 100 personas, 20 son fumadoras.
¿Cuál es la función de verosimilitud para la proporción de fumadores en la ciudad, \(\theta\)?
Llamemos \(Y\) al número de personas fumadoras en una muestra de tamaño \(n=100\). Suponiendo que la proporción de personas fumadoras en la ciudad es \(\theta\), \(Y\) sigue una distribución binomial con parámetros \(n=100\) y \(\theta\):
\[ Y|\theta \sim B(n=100, \theta), \]
por lo que la función de verosimilitud es
\[ L(\theta|Y=20) = P(Y=20|\theta) = \binom{100}{20} \theta^{20} (1-\theta)^{80}. \]library(bayesrules)
plot_binomial_likelihood(y = 20, n = 100) +
labs(x = "Proporción de fumadores (θ)", y = "Verosimilitud")
Una vez que tenemos la distribución a priori y la función de verosimilitud, podemos calcular la distribución a posteriori con la regla de Bayes:
\[ f(\theta|Y=20) = \frac{f(\theta)L(\theta|Y=20)}{P(Y=20)}. \]
Sin embargo, la constante normalizadora \(P(Y=20)\) es difícil de calcular, porque requiere acumular las probabilidades \(f(\theta)L(\theta|Y=20)\) para cualquier \(\theta\). Como \(\theta\) es una variable continua, esta suma se convierte en la integral
\[ P(Y=20) = \int_0^1 f(\theta)L(\theta|Y=20)\, d\theta. \]
Afortunadamente, gracias a la conjugación beta-binomial, sabemos que la distribución a posteriori de \(\theta\) dado \(Y=20\) es también una distribución beta, con parámetros \(\alpha=30+20=50\) y \(\beta=60+100-20=140\), es decir, \(\theta|Y=20 \sim \mbox{Beta}(50, 140)\).
plot_beta(alpha = 50, beta = 140) +
labs(x = "Proporción de fumadores (θ)", y = "Densidad")
Podemos dibujar las tres distribuciones, a priori, verosimilitud y a posteriori, en la misma gráfica con la función plot_beta_binomial() del paquete bayesrules.
plot_beta_binomial()
Dibuja la distribución beta a priori, la distribución binomial de la función de verosimilitud y la distribución beta a posteriori dado el número de éxitos observado en la muestra.
Parámetros:
alpha: parámetro \(\alpha\) de la distribución beta a priori.beta: parámetro \(\beta\) de la distribución beta a priori.y: número de éxitos observados en la muestra.n: tamaño de la muestra.plot_beta_binomial(alpha = 30, beta = 60, y = 20, n = 100) +
labs(x = "Proporción de fumadores (θ)", y = "Densidad / Verosimilitud")
Cuando no es posible determinar analíticamente la distribución a posteriori, como ocurre con modelos bayesianos con distribuciones no conjugadas, puede simularse mediante técnicas de muestreo. Las técnicas más habituales, por orden de aparición, son:
La técnica de muestreo de rejilla es la más sencilla y consiste en:
Ejemplo 4.9 Partiendo de la distribución a priori \(\theta \sim \mbox{Beta}(30, 60)\) y de la función de verosimilitud \(L(\theta|Y=20) = \binom{100}{20} \theta^{20} (1-\theta)^{80}\), el siguiente código simula la distribución a posteriori de \(\theta\) dado \(Y=20\) mediante muestreo de rejilla.
# Rejilla de 100 posibles valores de θ entre 0 y 1
rejilla <- seq(0, 1, length.out = 100)
# Función de verosimilitud para cada valor de θ
verosimilitud <- dbinom(20, size = 100, prob = rejilla)
# Probabilidad a priori para cada valor de θ
prior <- dbeta(rejilla, shape1 = 30, shape2 = 60)
# Plausibilidad a posteriori para cada valor de θ
posterior <- prior * verosimilitud
# Normalizamos la plausibilidad a posteriori para obtener una distribución de probabilidad
df_posterior <- tibble(θ = rejilla, posterior = posterior / sum(posterior))
# Gráfica de la distribución a posteriori
ggplot(df_posterior, aes(x = θ, y = posterior)) +
geom_line(color = myblue, linewidth = 1.5) +
labs(x = "Proporción de fumadores (θ)", y = "Probabilidad a posteriori (rejilla)") +
theme_minimal()
Para ver si la distribución simulada se parece a la distribución teórica \(\mbox{Beta}(50, 140)\), podemos extraer de ella una muestra aleatoria suficientemente grande y dibujar su histograma:
set.seed(123)
# Muestra aleatoria de tamaño 10000
muestra_posterior <- slice_sample(df_posterior, n = 10000, weight_by = posterior, replace = TRUE)
# Histograma de la muestra a posteriori
ggplot(muestra_posterior, aes(x = θ)) +
geom_histogram(aes(y = after_stat(density)), binwidth = 0.01, alpha = 0.7, fill = myblue, color = "white") +
stat_function(fun = dbeta, args = list(50, 140), color = myred, linewidth = 1) +
labs(x = "Proporción de fumadores (θ)", y = "Densidad") +
theme_minimal()
La distribución de frecuencias de la muestra simulada se ajusta bastante bien a la función de densidad de la distribución teórica \(\mbox{Beta}(50, 140)\).
Ejercicio 4.6 Partiendo de la distribución a posteriori anterior, \(\theta\sim \mbox{Beta}(50, 140)\), simula la distribución a posteriori si se toma otra muestra de tamaño 300 en la que aparecen 140 personas fumadoras.
# Rejilla de 100 posibles valores de θ entre 0 y 1
rejilla <- seq(0, 1, length.out = 100)
# Función de verosimilitud para cada valor de θ con la nueva muestra
verosimilitud <- dbinom(140, size = 300, prob = rejilla)
# Probabilidad a priori para cada valor de θ (distribución a posteriori anterior)
prior <- dbeta(rejilla, shape1 = 50, shape2 = 140)
# Plausibilidad a posteriori para cada valor de θ
posterior <- prior * verosimilitud
# Normalización de la plausibilidad a posteriori
df_posterior <- tibble(θ = rejilla, posterior = posterior / sum(posterior))
set.seed(123)
# Muestra aleatoria de tamaño 10000
muestra_posterior <- slice_sample(df_posterior, n = 10000, weight_by = posterior, replace = TRUE)
# Histograma de la muestra a posteriori
ggplot(muestra_posterior, aes(x = θ)) +
geom_histogram(aes(y = after_stat(density)), binwidth = 0.01, alpha = 0.7, fill = myblue, color = "white") +
stat_function(fun = dbeta, args = list(50 + 140, 140 + 300 - 140), color = myred, linewidth = 1) +
labs(x = "Proporción de fumadores (θ)", y = "Densidad") +
theme_minimal()
Cuando el parámetro de interés es una tasa \(\lambda\), como la tasa de nacimientos diarios en una ciudad, y los datos observados son el número de eventos \(X\) en un intervalo de tiempo dado, como el número de nacimientos en un día, el modelo bayesiano más utilizado es el modelo gamma-Poisson:
Ejercicio 4.7 Supongamos que queremos estimar la tasa de infartos diarios por cada millón de habitantes. Sabemos que en varios países similares la tasa media de infartos es de 5 por cada 1000 habitantes al año, y que la menor y la mayor tasas observadas fueron de 2 y 10 por cada 1000 habitantes al año, respectivamente.
¿Qué distribución gamma podríamos utilizar para modelizar la distribución a priori?
Puedes usar la aplicación distribuciones para probar distintas distribuciones gamma y encontrar la que mejor se ajuste a la información.
Una distribución \(\mbox{Gamma}(10, 2)\) tiene media \(10/2 = 5\) y un rango de valores compatible con la información disponible.
plot_gamma(shape = 10, rate = 2) +
labs(x = "Tasa de infartos (λ)", y = "Densidad")
Ejercicio 4.8 Tras realizar un nuevo sondeo, hemos observado que en una muestra de 5 días se han producido \(X_1=10, X_2=8, X_3=9, X_4=11, X_5=7\) infartos.
¿Cuál es la función de verosimilitud para la tasa de infartos \(\lambda\)?
Llamemos \(Y\) al número de infartos diarios. Suponiendo que la tasa de infartos diarios es \(\lambda\), el número de infartos de cada día sigue una distribución de Poisson de parámetro \(\lambda\), por lo que la función de verosimilitud para \(\lambda\) dados los cinco días observados es
\[ \begin{align*} L(\lambda|Y=(10,8,9,11,7)) = &P(Y_1=10|\lambda)P(Y_2=8|\lambda)P(Y_3=9|\lambda)\\ &P(Y_4=11|\lambda)P(Y_5=7|\lambda), \end{align*} \]
es decir,
\[ \begin{align*} L(\lambda|Y=(10,8,9,11,7)) &= \frac{e^{-\lambda}\lambda^{10}}{10!} \frac{e^{-\lambda}\lambda^{8}}{8!} \frac{e^{-\lambda}\lambda^{9}}{9!}\frac{e^{-\lambda}\lambda^{11}}{11!} \frac{e^{-\lambda}\lambda^{7}}{7!} \\ &= \frac{ e^{-5\lambda}\lambda^{45}}{10!\,8!\,9!\,11!\,7!}. \end{align*} \]plot_poisson_likelihood(y = c(10, 8, 9, 11, 7), lambda_upper_bound = 15) +
labs(x = "Tasa de infartos (λ)", y = "Verosimilitud")
Una vez que tenemos la distribución a priori y la función de verosimilitud, la distribución a posteriori se obtiene con la regla de Bayes:
\[ f(\lambda|Y=(10,8,9,11,7)) = \frac{f(\lambda)L(\lambda|Y=(10,8,9,11,7))}{P(Y=(10,8,9,11,7))}. \]
De nuevo, la constante normalizadora \(P(Y=(10,8,9,11,7))\) es difícil de calcular, porque requiere la integral
\[ P(Y=(10,8,9,11,7)) = \int_0^\infty f(\lambda)L(\lambda|Y=(10,8,9,11,7))\, d\lambda. \]
Pero, gracias a la conjugación gamma-Poisson, la distribución a posteriori de \(\lambda\) dada la muestra es también una distribución gamma con parámetros \(\alpha=10+45=55\) y \(\beta=2+5=7\), es decir, \(\lambda \sim \mbox{Gamma}(55, 7)\).
plot_gamma(shape = 55, rate = 7) +
labs(x = "Tasa de infartos (λ)", y = "Densidad")
Podemos dibujar las tres distribuciones, a priori, verosimilitud y a posteriori, en la misma gráfica con la función plot_gamma_poisson() del paquete bayesrules.
plot_gamma_poisson()
Dibuja la distribución gamma a priori, la distribución de Poisson de la función de verosimilitud y la distribución gamma a posteriori, dada una muestra de conteos.
Parámetros:
shape: parámetro \(\alpha\) de la distribución gamma a priori.rate: parámetro \(\beta\) de la distribución gamma a priori.sum_y: suma de los conteos observados en la muestra.n: tamaño de la muestra.plot_gamma_poisson(shape = 10, rate = 2, sum_y = 45, n = 5) +
labs(x = "Tasa de infartos (λ)", y = "Densidad / Verosimilitud")
Ejemplo 4.10 Partiendo de la distribución a priori \(\lambda \sim \mbox{Gamma}(10, 2)\) y de la función de verosimilitud \(L(\lambda|Y=(10, 8, 9, 11, 7))\), el siguiente código simula la distribución a posteriori de \(\lambda\) mediante muestreo de rejilla.
# Rejilla de 100 posibles valores de λ entre 1 y 15
rejilla <- seq(1, 15, length.out = 100)
# Función de verosimilitud para cada valor de λ
verosimilitud <- sapply(rejilla, function(l) prod(dpois(c(10, 8, 9, 11, 7), lambda = l)))
# Probabilidad a priori para cada valor de λ
prior <- dgamma(rejilla, shape = 10, rate = 2)
# Plausibilidad a posteriori para cada valor de λ
posterior <- prior * verosimilitud
# Normalizamos la plausibilidad a posteriori para obtener una distribución de probabilidad
df_posterior <- tibble(λ = rejilla, posterior = posterior / sum(posterior))
set.seed(123)
# Muestra aleatoria de tamaño 10000
muestra_posterior <- slice_sample(df_posterior, n = 10000, weight_by = posterior, replace = TRUE)
# Histograma de la muestra a posteriori
ggplot(muestra_posterior, aes(x = λ)) +
geom_histogram(aes(y = after_stat(density)), binwidth = 0.5, alpha = 0.7, fill = myblue, color = "white") +
stat_function(fun = dgamma, args = list(55, 7), color = myred, linewidth = 1) +
labs(x = "Tasa de infartos (λ)", y = "Densidad") +
theme_minimal()
Ejercicio 4.9 Partiendo de la distribución a posteriori anterior, \(\lambda\sim \mbox{Gamma}(55, 7)\), simula la distribución a posteriori si se toma otra muestra de tamaño 5 en la que se han producido \(X_1=5, X_2=7, X_3=4, X_4=6, X_5=3\) infartos.
# Rejilla de 100 posibles valores de λ entre 1 y 15
rejilla <- seq(1, 15, length.out = 100)
# Función de verosimilitud para cada valor de λ con la nueva muestra
verosimilitud <- sapply(rejilla, function(l) prod(dpois(c(5, 7, 4, 6, 3), lambda = l)))
# Probabilidad a priori para cada valor de λ (distribución a posteriori anterior)
prior <- dgamma(rejilla, shape = 55, rate = 7)
# Plausibilidad a posteriori para cada valor de λ
posterior <- prior * verosimilitud
# Normalización de la plausibilidad a posteriori
df_posterior <- tibble(λ = rejilla, posterior = posterior / sum(posterior))
set.seed(123)
muestra_posterior <- slice_sample(df_posterior, n = 10000, weight_by = posterior, replace = TRUE)
# Histograma de la muestra a posteriori
ggplot(muestra_posterior, aes(x = λ)) +
geom_histogram(aes(y = after_stat(density)), binwidth = 0.5, alpha = 0.7, fill = myblue, color = "white") +
stat_function(fun = dgamma, args = list(55 + 25, 7 + 5), color = myred, linewidth = 1) +
labs(x = "Tasa de infartos (λ)", y = "Densidad") +
theme_minimal()
Una de las propiedades más importantes de la regla de Bayes es que el orden en el que se obtienen los datos no afecta a la distribución a posteriori final:
\[ f(\theta|Y_1, Y_2) = f(\theta|Y_2, Y_1). \]
La demostración es sencilla. Si primero se observa \(Y_1\),
\[ f(\theta|Y_1) = \frac{f(\theta)L(\theta|Y_1)}{P(Y_1)}, \]
y después se observa \(Y_2\),
\[ \begin{align*} f(\theta|Y_1,Y_2) &= \frac{\frac{f(\theta)L(\theta|Y_1)}{P(Y_1)}L(\theta|Y_2)}{P(Y_2)} = \frac{f(\theta)L(\theta|Y_1)L(\theta|Y_2)}{P(Y_1)P(Y_2)} \\ &= \frac{f(\theta)L(\theta|Y_2)L(\theta|Y_1)}{P(Y_2)P(Y_1)} = f(\theta|Y_2,Y_1). \end{align*} \]
Cuando el parámetro de interés es una media \(\mu\), como el nivel medio de colesterol en sangre, y los datos \(X\) siguen una distribución normal con desviación típica conocida, el modelo bayesiano más utilizado es el modelo normal-normal:
Ejercicio 4.10 Supongamos que queremos estimar el nivel medio de colesterol en sangre en una población. Sabemos que en poblaciones similares el nivel medio de colesterol está en un rango que va de 191 a 209 mg/dL.
¿Qué distribución normal podríamos utilizar para modelizar la distribución a priori?
Puedes usar la aplicación distribuciones para probar distintas distribuciones normales y encontrar la que mejor se ajuste a la información.
Una distribución \(N(200, 3^2)\) tiene media 200 y desviación típica 3, de modo que aproximadamente el 95 % de su masa cae entre 194 y 206 mg/dL, dentro del rango de la información disponible.
library(bayesrules)
plot_normal(mean = 200, sd = 3) +
labs(x = "Nivel de colesterol (mg/dL)", y = "Densidad")
Ejercicio 4.11 Tras realizar un nuevo sondeo, hemos observado que en una muestra de 10 personas el nivel de colesterol en sangre es 230, 226, 245, 228, 232, 235, 232, 244, 251 y 247 mg/dL. Suponiendo que el nivel de colesterol sigue una distribución normal con desviación típica conocida de 40 mg/dL, algo razonable en una variable biológica continua, ¿cuál es la función de verosimilitud para el nivel medio de colesterol \(\mu\)?
Llamemos \(Y\) al nivel de colesterol en sangre. Si \(Y\sim N(\mu, 40^2)\), la función de verosimilitud para \(\mu\) dados los valores observados \(y=(y_1,\ldots,y_n)\) es
\[ L(\mu|y) = \frac{1}{(2\pi \sigma^2)^{n/2}} \exp\left(-\frac{1}{2\sigma^2} \sum_{i=1}^n (y_i - \mu)^2\right), \qquad \sigma^2 = 40^2. \]muestra <- c(230, 226, 245, 228, 232, 235, 232, 244, 251, 247)
plot_normal_likelihood(y = muestra, sigma = 40) +
labs(x = "Nivel medio de colesterol (μ)", y = "Verosimilitud")
Una vez definidas la verosimilitud y la distribución a priori, la distribución a posteriori se obtiene con la regla de Bayes. La función plot_normal_normal() del paquete bayesrules dibuja las tres distribuciones en la misma gráfica:
plot_normal_normal()
Dibuja la distribución normal a priori, la distribución normal de la función de verosimilitud y la distribución normal a posteriori, dada una muestra de datos.
Parámetros:
mean: media de la distribución normal a priori.sd: desviación típica de la distribución normal a priori.y_bar: media muestral.n: tamaño de la muestra.sigma: desviación típica conocida de la distribución normal de la función de verosimilitud.media <- mean(c(230, 226, 245, 228, 232, 235, 232, 244, 251, 247))
plot_normal_normal(mean = 200, sd = 3, y_bar = media, n = 10, sigma = 40) +
labs(x = "Nivel medio de colesterol (μ)", y = "Densidad / Verosimilitud")
El muestreo de rejilla es sencillo, pero no es factible para modelos con muchos parámetros, porque el número de puntos de la rejilla crece exponencialmente con el número de parámetros. Con 100 puntos por parámetro:
| Número de parámetros | Número total de puntos en la rejilla |
|---|---|
| 1 | 100 |
| 2 | 10.000 |
| 3 | 1.000.000 |
| 4 | 100.000.000 |
| 5 | 10.000.000.000 |
Por esta razón, para modelos con muchos parámetros se utilizan técnicas de muestreo más eficientes, como MCMC o HMC.
Los métodos de Monte Carlo generan muestras representativas de la distribución a posteriori con un coste computacional mucho menor que el muestreo de rejilla. Los dos tipos principales son:
Stan es un lenguaje de programación específico para el modelado y la inferencia bayesiana que implementa algoritmos de muestreo eficientes, como HMC.
Un programa en Stan se compone de varios bloques declarativos:
data: declara las variables observadas en los datos, sus dimensiones, tipos y rango posible de valores.parameters: declara los parámetros del modelo que queremos estimar, sus dimensiones, tipos y rango posible de valores.model: define la estructura del modelo bayesiano, especificando la distribución a priori de los parámetros y la función de verosimilitud de los datos dados los parámetros.generated quantities: permite definir variables adicionales que se calculan a partir de los parámetros y los datos, como predicciones o medidas de ajuste del modelo.Ejemplo 4.11 A continuación se muestra un programa en Stan para el modelo beta-binomial de la proporción de fumadores, con \(Y=20\) fumadoras en una muestra de 100 personas.
// modelo-beta-binomial.stan
// Datos
data {
int<lower = 0, upper = 100> Y;
}
// Parámetros
parameters {
real<lower = 0, upper = 1> pi;
}
// Modelo
model {
Y ~ binomial(100, pi);
pi ~ beta(30, 60);
}
// Posterior predictivo: simular nuevas observaciones
generated quantities {
int y_rep = binomial_rng(100, pi);
}rstanEl paquete rstan es la interfaz de R para Stan. Permite compilar y ejecutar modelos escritos en Stan desde R, así como analizar los resultados. Los programas se escriben como una cadena de caracteres y se compilan y ejecutan con la función stan().
stan()
Compila un programa de Stan y simula el muestreo de la distribución a posteriori dados los datos.
model_code: cadena de caracteres con el código del programa de Stan.data: lista con los datos que se pasan al programa.chains: número de cadenas de Markov que se simulan.iter: número de iteraciones de cada cadena.seed: semilla para la generación de números aleatorios.Ejemplo 4.12 El siguiente código define un modelo gamma-Poisson para la tasa de infartos diarios, con una distribución a priori \(\lambda \sim \mbox{Gamma}(10, 2)\) y una función de verosimilitud de Poisson. Lo compila y simula el muestreo de la distribución a posteriori.
library(rstan)
# Código del programa de Stan
modelo <- "
data {
int<lower=0> n; // número de días
int<lower=0> y[n]; // número de infartos diarios
}
parameters {
real<lower=0> lambda; // tasa de infartos
}
model {
y ~ poisson(lambda); // función de verosimilitud
lambda ~ gamma(10, 2); // distribución a priori
}
"
tasa_infartos_post <- stan(model_code = modelo,
data = list(n = 5, y = c(10, 8, 9, 11, 7)),
chains = 4, iter = 5000*2, seed = 123)El resultado es una muestra de 4 cadenas de Markov de 5000 valores cada una, es decir, un total de 20.000 muestras de la distribución a posteriori de \(\lambda\) dado el número de infartos diarios observado.
as.data.frame(tasa_infartos_post) |> head() lambda lp__
1 8.450915 58.22870
2 6.963386 57.99292
3 8.179375 58.33324
4 7.334871 58.25109
5 7.943305 58.37498
6 8.568749 58.16545Podemos visualizar las trazas de las cuatro cadenas con la función mcmc_trace del paquete bayesplot.
library(bayesplot)
mcmc_trace(tasa_infartos_post, pars = "lambda", size = 0.2) +
labs(x = "Iteración", y = "Valor de λ") +
theme_gray()
Si hacemos zoom sobre una de estas cadenas, podemos observar cómo recorre los distintos valores de la distribución a posteriori.
mcmc_trace(as.array(tasa_infartos_post)[, 1, , drop = FALSE], pars = "lambda", window = c(1, 100)) +
labs(x = "Iteración", y = "Valor de λ") +
theme_gray()
Podemos visualizar la distribución a posteriori de \(\lambda\) con la función mcmc_dens del paquete bayesplot, y compararla con la distribución teórica \(\mbox{Gamma}(55, 7)\):
mcmc_dens(tasa_infartos_post, pars = "lambda") +
stat_function(fun = dgamma, args = list(55, 7), color = myred) +
labs(x = "Tasa de infartos (λ)", y = "Densidad") +
theme_gray()
Ejercicio 4.12 Define en Stan un modelo normal-normal para el nivel medio de colesterol en sangre, con una distribución a priori \(N(200, 3^2)\) y una función de verosimilitud normal con desviación típica conocida de 40 mg/dL. Compílalo y simula el muestreo de la distribución a posteriori dada la muestra de 10 personas del ejercicio anterior.
# Código del programa de Stan
modelo <- "
data {
int<lower=0> n; // número de personas
vector[n] y; // niveles de colesterol
}
parameters {
real mu; // nivel medio de colesterol
}
model {
y ~ normal(mu, 40); // función de verosimilitud
mu ~ normal(200, 3); // distribución a priori
}
"
colesterol_post <- stan(model_code = modelo,
data = list(n = 10, y = c(230, 226, 245, 228, 232, 235, 232, 244, 251, 247)),
chains = 4, iter = 5000*2, seed = 123)mcmc_dens(colesterol_post, pars = "mu") +
stat_function(fun = dnorm, args = list(201.97, sqrt(8.52)), color = myred) +
labs(x = "Nivel medio de colesterol (μ)", y = "Densidad") +
theme_gray()
El principal problema de las muestras generadas con MCMC o HMC es que los valores de una cadena no son independientes: cada valor depende del anterior. Aunque el algoritmo busca que la muestra se ajuste a la distribución a posteriori, a veces puede fallar y generar muestras poco representativas.
Para analizar la calidad de las muestras se utilizan, entre otras, las siguientes técnicas:
Para analizar la distribución de frecuencias de cada cadena de Markov podemos usar la función mcmc_dens_overlay del paquete bayesplot.
mcmc_dens_overlay()
Dibuja las distribuciones de frecuencias de cada cadena de Markov generada.
Parámetros:
fit: objeto de clase stanfit con el resultado del muestreo de la distribución a posteriori.pars: vector de caracteres con los nombres de los parámetros que se analizan.Ejemplo 4.13
mcmc_dens_overlay(tasa_infartos_post, pars = "lambda") +
labs(x = "Tasa de infartos (λ)", y = "Densidad") +
theme_gray()
Cada valor de la cadena de Markov depende del anterior, lo que genera una cadena de dependencias \(\theta_1 \to \theta_2 \to \theta_3 \to \ldots \to \theta_n\). Es de esperar que la dependencia entre \(\theta_{i}\) y \(\theta_{i+k}\) disminuya a medida que \(k\) aumenta. Esta dependencia entre los valores de la cadena se conoce como autocorrelación.
Definición 4.3 (Autocorrelación de una cadena de Markov) Dada una cadena de Markov \(\theta_i\) con \(i=1, \ldots, n\), la autocorrelación de orden \(k\) se define como
\[ \rho_k = \frac{\text{Cov}(\theta_i, \theta_{i+k})}{\text{Var}(\theta_i)}, \]
donde \(\text{Cov}(\theta_i, \theta_{i+k})\) es la covarianza entre el valor \(\theta_i\) y el valor \(\theta_{i+k}\), y \(\text{Var}(\theta_i)\) es la varianza de los valores \(\theta_i\).
Para representar la autocorrelación de las muestras generadas podemos usar la función mcmc_acf del paquete bayesplot.
mcmc_acf()
Dibuja la función de autocorrelación de las muestras generadas por MCMC o HMC.
Parámetros:
fit: objeto de clase stanfit con el resultado del muestreo de la distribución a posteriori.pars: vector de caracteres con los nombres de los parámetros que se analizan.lags: número de retardos para calcular la autocorrelación.Ejemplo 4.14
mcmc_acf(tasa_infartos_post, pars = "lambda", lags = 10) +
labs(x = "Orden o retardo (k)", y = "Autocorrelación (ρ_k)") +
theme_gray()
La eficiencia de una cadena de Markov se puede evaluar calculando el número efectivo de muestras independientes (ESS), es decir, el número de muestras independientes a las que equivale la cadena autocorrelada.
Podemos analizar la eficiencia de las muestras con la función neff_ratio del paquete bayesplot.
neff_ratio()
Calcula la proporción entre el número efectivo de muestras independientes (ESS) y el número total de muestras generadas por MCMC o HMC.
Parámetros:
fit: objeto de clase stanfit con el resultado del muestreo de la distribución a posteriori.pars: vector de caracteres con los nombres de los parámetros que se analizan.Ejemplo 4.15
neff_ratio(tasa_infartos_post, pars = "lambda")[1] 0.3731105Cuanto más cercano a 1 sea este cociente, más eficiente es el muestreo.
La consistencia de las muestras se puede evaluar comparando la variabilidad de los valores muestreados en todas las cadenas combinadas con la variabilidad dentro de cada cadena individual.
Definición 4.4 (Métrica \(\hat{R}\)) Dada una colección de cadenas de Markov generadas por un mismo proceso (MCMC o HMC), la métrica \(\hat{R}\) es la raíz cuadrada de la razón entre la variabilidad total de los valores muestreados en todas las cadenas combinadas y la variabilidad dentro de cada cadena individual:
\[ \hat{R} = \sqrt{\frac{\text{Var}_{total}}{\text{Var}_{individual}}}. \]
Si la variabilidad entre cadenas es similar a la variabilidad dentro de cada cadena, \(\hat{R} \approx 1\), lo que indica que las cadenas han convergido a la misma distribución a posteriori y las muestras son consistentes.
Podemos calcular \(\hat{R}\) con la función rhat del paquete bayesplot.
rhat()
Calcula la métrica \(\hat{R}\) para evaluar la consistencia de las muestras generadas por MCMC o HMC.
Parámetros:
fit: objeto de clase stanfit con el resultado del muestreo de la distribución a posteriori.pars: vector de caracteres con los nombres de los parámetros que se analizan.Ejemplo 4.16
rhat(tasa_infartos_post, pars = "lambda")[1] 1.000579Si el valor es prácticamente 1, las cuatro cadenas de Markov generadas son consistentes entre sí.
Ejercicio 4.13 Analiza la traza, la autocorrelación, la eficiencia y la consistencia de las cadenas de Markov generadas con Stan para la distribución a posteriori del nivel medio de colesterol en sangre del ejercicio anterior.
# Análisis de la traza
mcmc_trace(colesterol_post, pars = "mu", size = 0.2) +
labs(x = "Iteración", y = "Valor de μ")
# Análisis de la autocorrelación
mcmc_acf(colesterol_post, pars = "mu", lags = 10) +
labs(x = "Orden o retardo (k)", y = "Autocorrelación (ρ_k)")
# Eficiencia
neff_ratio(colesterol_post, pars = "mu")[1] 0.3672184# Consistencia
rhat(colesterol_post, pars = "mu")[1] 1.000358Una de las críticas más comunes a la inferencia bayesiana es que los resultados pueden depender de la elección de la distribución a priori, lo que puede introducir un sesgo subjetivo en el análisis.
Los defensores de la estadística bayesiana argumentan que la elección de la distribución a priori debe basarse en el conocimiento previo disponible sobre el fenómeno estudiado, y no en una preferencia personal o arbitraria. En cualquier caso, es importante realizar un análisis de sensibilidad: evaluar cómo cambian los resultados a posteriori al cambiar la distribución a priori. Así se identifica si los resultados son robustos o si dependen en gran medida de la elección.
Una forma de reducir la influencia de la distribución a priori es usar distribuciones poco informativas, que asignan una probabilidad similar a un amplio rango de valores del parámetro. Como ejemplo, consideremos la estimación del nivel medio de colesterol con una distribución a priori uniforme continua \(U(150, 300)\).
Ejemplo 4.17 Veamos cómo cambia la distribución a posteriori del nivel medio de colesterol con la muestra de 10 personas, usando una distribución a priori uniforme continua.
# Código del programa de Stan
modelo <- "
data {
int<lower=0> n; // número de personas
vector[n] y; // niveles de colesterol
}
parameters {
real<lower=150, upper=300> mu; // nivel medio de colesterol
}
model {
y ~ normal(mu, 40); // función de verosimilitud
mu ~ uniform(150, 300); // distribución a priori
}
"
posterior_sim <- stan(model_code = modelo,
data = list(n = 10, y = c(230, 226, 245, 228, 232, 235, 232, 244, 251, 247)),
chains = 4, iter = 5000*2, seed = 123)mcmc_dens(posterior_sim, pars = "mu") +
stat_function(fun = dnorm, args = list(201.97, sqrt(8.52)), color = myred) +
labs(x = "Nivel medio de colesterol (μ)", y = "Densidad") +
theme_gray()
Con una muestra más grande, de 100 personas, simulada a partir de la misma distribución, la distribución a posteriori es prácticamente la misma con cualquier prior razonable:
set.seed(123)
posterior_sim <- stan(model_code = modelo,
data = list(n = 10^2, y = rnorm(10^2, mean = 202, sd = 40)),
chains = 4, iter = 5000*2, seed = 123)mcmc_dens(posterior_sim, pars = "mu") +
labs(x = "Nivel medio de colesterol (μ)", y = "Densidad") +
theme_gray()
Para muestras pequeñas puede haber bastante discrepancia entre las distribuciones a posteriori generadas con distintas distribuciones a priori, pero esa diferencia se disipa a medida que tenemos más datos.
Cuando la muestra es pequeña, la elección de la distribución a priori es determinante. Según el grado de información que incorporan, se distinguen tres tipos:
Una distribución a priori no informativa asigna la misma plausibilidad a todo el rango de valores del parámetro. El caso más extremo es el modelo uniforme impropio \(f(\theta) \propto 1\), que no es una distribución de probabilidad porque no integra 1, pero que produce distribuciones a posteriori adecuadas cuando hay datos suficientes.
¿Cuándo usarla? Cuando no tenemos información previa y los datos son abundantes.
¿Cuándo no usarla? No debería usarse con muestras pequeñas, porque la posterior puede ser inestable o producir estimaciones absurdas.
Ejemplo 4.18 Para estimar la tensión arterial sistólica media de una población se utiliza un modelo uniforme \(f(\mu) \propto 1\), que asigna la misma plausibilidad a cualquier valor de \(\mu>0\).
Una distribución a priori poco informativa asigna una probabilidad similar a un amplio rango de valores del parámetro, pero descarta los valores que no tienen sentido. Suele utilizarse una distribución uniforme continua con el rango restringido a los valores plausibles, o una distribución normal con desviación típica grande para reflejar la incertidumbre sobre el parámetro.
¿Cuándo usarla? Cuando no tenemos información previa pero queremos evitar que la posterior sea inestable con muestras pequeñas. Es la opción más habitual cuando no hay información sólida sobre la incertidumbre del parámetro.
¿Cuándo no usarla? No debería usarse cuando tenemos información previa sólida, porque no aprovecha esa información para mejorar la precisión de las estimaciones, especialmente con pocos datos.
Ejemplo 4.19 Para estimar la tensión arterial sistólica media de una población se utiliza un modelo uniforme \(\mu \sim U(70, 200)\), que asigna la misma plausibilidad a cualquier valor de \(\mu\) entre 70 y 200 mmHg, y descarta los valores menores de 70 o mayores de 200 mmHg, que no tienen sentido en este contexto.
Una distribución a priori informativa refleja un conocimiento previo sólido sobre el parámetro de interés (resultados de estudios anteriores, valores de referencia clínica u opinión de expertos). Puede mejorar la precisión de las estimaciones cuando hay poca información en los datos. Suelen elegirse distribuciones conjugadas con la distribución de los datos.
¿Cuándo usarla? Cuando existe evidencia previa sólida (metaanálisis, registros históricos, valores fisiológicos bien establecidos). Siempre debe documentarse y justificarse.
¿Cuándo no usarla? No debería usarse cuando no existe información previa sólida, porque puede introducir un sesgo subjetivo y producir estimaciones poco fiables, especialmente con pocos datos.
Ejemplo 4.20 Para estimar la tensión arterial sistólica media de una población se utiliza un modelo normal \(\mu \sim N(120, 15^2)\), que refleja el conocimiento previo sobre la tensión arterial sistólica media en poblaciones similares, con una media de 120 mmHg y una desviación típica de 15 mmHg para reflejar la incertidumbre sobre ese valor.
| Situación | Prior recomendado |
|---|---|
| Media poblacional | normal(valor_referencia, incertidumbre) |
| Desviación típica | half_t(3, 0, 1) o exponential(1) |
| Proporción sin información previa | beta(1, 1) o beta(0.5, 0.5) |
| Proporción con información previa | beta() calibrada con la media y la varianza esperadas |
| Coeficientes de regresión (datos estandarizados) | normal(0, 1) |
| Coeficientes de regresión (datos sin estandarizar) | normal(0, 2.5), recomendación de Gelman |
| Muchos predictores, selección de variables | horseshoe(df=1) |
| Regresión logística, coeficientes | normal(0, 2.5) o student_t(7, 0, 2.5) |
| Parámetro | Distribución prior | Parámetros habituales | Conjugada con | Aplicación en epidemiología |
|---|---|---|---|---|
| Prevalencia, proporción | Beta(\(\alpha, \beta\)) | Beta(1,1) uniforme · Beta(0,5;0,5) Jeffreys · Beta(2,5) informativa | Binomial / Bernoulli | Prevalencia de enfermedad, tasa de ataque, sensibilidad y especificidad de un test |
| Vector de proporciones | Dirichlet(\(\alpha_1,\ldots,\alpha_k\)) | Dir(1,…,1) uniforme · Dir(0,5,…,0,5) Jeffreys | Multinomial / Categórica | Distribución de serotipos, variantes, grupos de riesgo |
| Tasa de incidencia, tasa de mortalidad | Gamma(\(\alpha, \beta\)) | Gamma(0,01; 0,01) difusa · Gamma(2; 0,5) débil | Poisson / Exponencial | Incidencia de casos por 100.000 habitantes, mortalidad, tasa de hospitalización |
| Número reproductivo \(R_0\), \(R(t)\) | Gamma(\(\alpha, \beta\)) | Gamma(1; 0,2) · Gamma(2; 0,5) | Poisson (ecuación de renovación) | Transmisibilidad de enfermedades infecciosas, estimación de \(R(t)\) |
| Media poblacional (tensión arterial, colesterol, glucemia) | Normal(\(\mu_0, \sigma_0^2\)) | N(valor de referencia, \(\sigma\) amplio) | Normal (\(\sigma^2\) conocida) | Estimación de la media de biomarcadores clínicos |
| Media y varianza poblacional | Normal-inversa-gamma(\(\mu_0,\lambda,\alpha,\beta\)) | Hiperparámetros basados en valores de referencia | Normal completa | Biomarcadores con media y varianza desconocidas |
| Log-odds, coeficiente de regresión logística | Normal(\(0, \sigma^2\)) | N(0; 2,5²) de Gelman · N(0, 1) estandarizado | — (MCMC) | Factores de riesgo en estudios de casos y controles, OR ajustado |
| Log-hazard ratio, coeficiente de supervivencia | Normal(\(0, \sigma^2\)) | N(0; 2,5²) | — (MCMC) | Análisis de supervivencia: tiempo hasta el evento, censura |
| Desviación típica, varianza residual | half-normal(\(0, \sigma\)) · half-\(t(\nu,0,\sigma)\) | half-N(0,1) · half-\(t(3,0,1)\) de Gelman | — (MCMC) | Variabilidad residual en regresión, modelos mixtos |
| Varianza de efectos aleatorios | half-\(t(\nu, 0, \sigma)\) · gamma inversa(\(\alpha,\beta\)) | half-\(t(3,0,1)\) recomendado · IG(0,01; 0,01) clásica | — (MCMC) | Heterogeneidad entre centros, hospitales o regiones en modelos jerárquicos |
| Correlación entre efectos aleatorios | LKJ(\(\eta\)) | LKJ(1) uniforme · LKJ(2) concentrada | — (MCMC) | Correlación entre interceptos y pendientes en modelos mixtos |
| Eficacia vacunal (escala log-OR) | Normal(\(0, \sigma^2\)) · Beta (escala de probabilidad) | N(\(\log(0{,}5)\), 0,5²) informativa | — (MCMC) | Ensayos clínicos de vacunas, estudios de efectividad vacunal |
| Sensibilidad y especificidad de un test | Beta(\(\alpha, \beta\)) | Beta(18, 2) para Se ≈ 90 % · Beta(38, 2) para Sp ≈ 95 % | Binomial | Evaluación de pruebas diagnósticas, modelos de clasificación errónea |
| Coeficientes con muchas variables | horseshoe(df) · Normal(0,1) | horseshoe(df=1) · N(0,1) estandarizado | — (MCMC) | Selección de factores de riesgo, estudios de exposición múltiple |
| Parámetro de forma de Weibull (supervivencia) | Gamma(\(\alpha, \beta\)) · half-normal | Gamma(1,1) · half-N(0,1) | — (MCMC) | Modelos paramétricos de supervivencia, tiempo hasta recaída o muerte |