Introducción

El modelo Normal para variables continuas \(y_i \in \mathbb{R}\), con \(i = 1, \dots, n\), se define como:
\[ \begin{aligned} y_i \mid \theta, \sigma^2 &\overset{\text{iid}}{\sim} \textsf{N}(\theta, \sigma^2) \\ (\theta, \sigma^2) &\sim p(\theta, \sigma^2) \end{aligned} \] donde \(\boldsymbol{\theta} = (\theta, \sigma^2) \in \Theta = \mathbb{R} \times \mathbb{R}^+\).

Considere especificar el estado de información previo sobre \(\theta\) de manera independiente de \(\sigma^2\), de modo que \(p(\theta,\sigma^2)=p(\theta)\,p(\sigma^2)\), donde \[ \theta\sim\textsf{N}(\mu_0,\tau^2_0)\qquad\text{y}\qquad\sigma^2\sim\textsf{GI}\left(\frac{\nu_0}{2},\frac{\nu_0\sigma^2_0}{2}\right), \]
siendo \(\mu_0\), \(\tau^2_0\), \(\nu_0\) y \(\sigma^2_0\) los hiperparámetros del modelo.

Esta formulación es flexible ya que no impone una restricción de dependencia previa entre \(\theta\) y \(\sigma^2\).

En el caso de una previa conjugada, donde \(\tau_0^2\) es proporcional a \(\sigma^2\), la distribución posterior de \(\sigma^2\) sigue una Gamma Inversa. Sin embargo, en este caso, la distribución posterior de \(\sigma^2\) no corresponde a una distribución estándar.

Muetreador de Gibbs

El muestreador de Gibbs (Gibbs sampler) es un algoritmo que permite generar muestras correlacionadas de la distribución posterior utilizando las distribuciones condicionales completas de los parámetros.

Dado un estado actual de los parámetros del modelo \(\boldsymbol{\theta}^{(b)} = (\theta^{(b)},\sigma^{2\,(b)})\), se genera un nuevo estado \(\boldsymbol{\theta}^{(b+1)} = (\theta^{(b+1)},\sigma^{2\,(b+1)})\) mediante el siguiente procedimiento iterativo:

Muestrear \(\theta^{(b+1)}\sim p(\theta\mid\sigma^{2\,(b)},\boldsymbol{y})\).
Muestrear \(\sigma^{2\,(b+1)}\sim p(\sigma^2\mid\theta^{(b+1)},\boldsymbol{y})\).
Actualizar \(\boldsymbol{\theta}^{(b+1)} = (\theta^{(b+1)},\sigma^{2\,(b+1)})\).
Repetir los pasos 1 a 3 hasta alcanzar la convergencia.

Este procedimiento genera una secuencia dependiente de muestras \(\boldsymbol{\theta}^{(1)}, \ldots, \boldsymbol{\theta}^{(B)}\), que provienen de la distribución posterior \(p(\theta, \sigma^2 \mid \boldsymbol{y})\).

El muestreador de Gibbs produce una cadena de Markov (Markov chain) \(\boldsymbol{\theta}^{(1)},\dots,\boldsymbol{\theta}^{(B)}\), ya que cada nuevo estado \(\boldsymbol{\theta}^{(b+1)}\) depende exclusivamente del estado anterior \(\boldsymbol{\theta}^{(b)}\).

La distribución estacionaria de la cadena \(\boldsymbol{\theta}^{(1)},\dots,\boldsymbol{\theta}^{(B)}\) es la distribución posterior \(p(\boldsymbol{\theta} \mid \boldsymbol{y})\).

Referencias:

Smith, A. F. M., & Roberts, G. O. (1993). Bayesian computation via the Gibbs sampler and related Markov chain Monte Carlo methods. Journal of the Royal Statistical Society: Series B (Methodological), 55(1), 3–23.

Distribuciones condicionales completas

Bajo el modelo Normal con distribución previa semi-conjugada, se tiene que:

(Ejercicio.) La distribución condicional completa de \(\theta\) es \(\theta\mid\sigma^2,\boldsymbol{y}\sim\textsf{N}(\mu_n,\tau^2_n)\), donde \[ \mu_n = \frac{\frac{1}{\tau^2_0}\mu_0 + \frac{n}{\sigma^2}\bar{y}}{\frac{1}{\tau^2_0} + \frac{n}{\sigma^2}} \qquad\text{y}\qquad\tau^2_n=\frac{1}{\frac{1}{\tau^2_0} + \frac{n}{\sigma^2}}\,. \]
(Ejercicio.) La distribución condicional completa de \(\sigma^2\) es \(\sigma^2\mid\theta,\boldsymbol{y}\sim\textsf{GI}\left(\tfrac{\nu_n}{2},\tfrac{\nu_n\,\sigma^2_n}{2}\right)\), donde \[ \nu_n = \nu_0+n\qquad\text{y}\qquad \nu_n\sigma^2_n = \nu_0\sigma^2_0 + ns^2_\theta\,, \] con \(s^2_\theta = \tfrac{1}{n}\sum_{i=1}^n (y_i-\theta)^2\), lo que corresponde a un estimador insesgado de \(\sigma^2\) si \(\theta\) fuera conocido.

Como punto de partida, basta con inicializar a \(\sigma^2\). Usualmente, este valor se genera a partir de la distribución previa, es decir, \(\sigma^{2\,(0)}\sim\textsf{IG}\left(\tfrac{\nu_0}{2},\tfrac{\nu_0\sigma^2_0}{2}\right)\).

Ejemplo: Puntajes de Matemáticas

La base de datos contiene la información de una muestra aleatoria simple de los estudiantes que presentaron la Prueba Saber 11 en 2023-2.

La prueba de matemáticas está diseñada con una escala de 0 a 100 (sin decimales), un puntaje promedio de 50 y una desviación estándar de 10 puntos.

¿Existe suficiente evidencia para concluir que el puntaje promedio de Matemáticas de Bogotá es significativamente superior al puntaje promedio preestablecido por la prueba?

Los datos son públicos y están disponibles en este enlace.

Tratamiento de datos

# Cargar datos
dat <- read.csv("SB11_20232_muestra.txt", sep = ";", stringsAsFactors = FALSE)

# Dimensiones
dim(dat)

## [1] 5511   83

# Distribución de frecuencias de la variable ESTU_DEPTO_RESIDE
table(dat$ESTU_DEPTO_RESIDE)

## 
##        AMAZONAS       ANTIOQUIA          ARAUCA       ATLANTICO          BOGOTÁ 
##               8             758              34             342             776 
##         BOLIVAR          BOYACA          CALDAS         CAQUETA        CASANARE 
##             275             173              78              42              53 
##           CAUCA           CESAR           CHOCO         CORDOBA    CUNDINAMARCA 
##             126             149              55             218             389 
##         GUAINIA        GUAVIARE           HUILA      LA GUAJIRA       MAGDALENA 
##               4              15             126             100             200 
##            META          NARIÑO NORTE SANTANDER        PUTUMAYO         QUINDIO 
##             122             145             147              38              45 
##       RISARALDA      SAN ANDRES       SANTANDER           SUCRE          TOLIMA 
##             125               9             274             107             170 
##           VALLE          VAUPES         VICHADA 
##             401               3               4

# Datos Bogotá
y <- dat[dat$ESTU_COD_RESIDE_DEPTO == 11, "PUNT_MATEMATICAS"]

# Tamaño de la muestra
(n <- length(y))

## [1] 776

# Estadísticos suficientes
(yb <- mean(y))

## [1] 55.08247

(s2 <- var(y))

## [1] 143.8487

Muestreador de Gibbs

# Hiperparámetros
mu0 <- 50  
t20 <- 10^2  # Regla empírica: P( |\theta - \mu_0| < 3\tau_0 ) = 99.7%
s20 <- 10^2 
nu0 <- 1

# Número de muestras 
B <- 10000

# Matriz para almacenar las muestras
THETA <- matrix(data = NA, nrow = B, ncol = 2)

# Algoritmo: Muestreador de Gibbs

# Inicialización
set.seed(123)
sig2 <- 1 / rgamma(1, shape = nu0 / 2, rate = nu0 * s20 / 2)  # Valor inicial de sigma^2

# Iteraciones del muestreador de Gibbs
set.seed(123)
for (b in 1:B) {
  # Actualización de theta
  t2n   <- 1 / (1 / t20 + n / sig2)
  mun   <- t2n * (mu0 / t20 + n * yb / sig2)
  theta <- rnorm(1, mean = mun, sd = sqrt(t2n))
  
  # Actualización de sigma^2
  nun  <- nu0 + n 
  s2n  <- (nu0 * s20 + sum((y - theta)^2)) / nun
  sig2 <- 1 / rgamma(1, shape = nun / 2, rate = nun * s2n / 2)
  
  # Almacenar muestras
  THETA[b, ] <- c(theta, sig2)
  
  # Mostrar progreso cada 10% de las iteraciones
  if (b %% floor(B / 10) == 0) 
    cat(sprintf("%.0f%% completado...\n", 100 * b / B))
}

## 10% completado...
## 20% completado...
## 30% completado...
## 40% completado...
## 50% completado...
## 60% completado...
## 70% completado...
## 80% completado...
## 90% completado...
## 100% completado...

# Establecer THETA como data frame
THETA <- as.data.frame(THETA)
colnames(THETA) <- c("theta", "sigma2")

Cadenas

# Cadenas
par(mfrow = c(2,1), mar = c(2.75, 3, 1.5, 0.5), mgp = c(1.7, 0.7, 0))

# Gráfico de la cadena para theta
plot(THETA$theta[-c(1:10)], type = "p", pch = 16, col = adjustcolor(1, 0.3), cex = 0.8, 
     xlab = "Iteración", ylab = expression(theta), 
     main = expression("Cadena de " ~ theta))

# Gráfico de la cadena para sigma^2
plot(THETA$sigma2[-c(1:10)], type = "p", pch = 16, col = adjustcolor(1, 0.3), cex = 0.8, 
     xlab = "Iteración", ylab = expression(sigma^2), 
     main = expression("Cadena de " ~ sigma^2))

Distribución posterior

Inferencia sobre la media y la varianza.
	Estimación	CV	L. Inf. 95%	L. Sup. 95%
Media	55.069	0.008	54.230	55.900
Varianza	144.200	0.051	130.521	159.183

Distriución predictiva

# Distribución predictiva posterior
set.seed(1234)
y_new <- rnorm(n = B, mean = THETA$theta, sd = sqrt(THETA$sigma2))

Inferencia sobre una predicción.
	Estimación	CV	L. Inf. 95%	L. Sup. 95%
y pred	55.137	0.215	31.858	78.424

Chequeo del modelo

# Estadísticos observados
t_obs <- c(mean(y), sd(y))

# Inicializar matriz para almacenar estadísticas de prueba
t_mc <- matrix(NA, nrow = B, ncol = 2)

# Distribución predictiva posterior
set.seed(1234)
for (i in 1:B) {
  # Datos simulados
  y_rep <- rnorm(n = n, mean = THETA$theta[i], sd = sqrt(THETA$sigma2[i]))

  # Estadísticos de prueba
  t_mc[i, ] <- c(mean(y_rep), sd(y_rep))
}

# ppp
ppp_media <- round(mean(t_mc[ , 1] < t_obs[1]), 3)
ppp_sd    <- round(mean(t_mc[ , 2] < t_obs[2]), 3)

Inferencia: Bogotá

Inferencia sobre el puntaje promedio de Matemáticas en Bogotá.
	Estimación	CV	L. Inf. 95%	L. Sup. 95%
Bayesiana	55.069	0.008	54.230	55.900
Frec. Asintótica	55.082	0.008	54.239	55.926
Frec. Bootstrap Par.	55.079	0.008	54.232	55.922
Frec. Bootstrap No Par.	55.080	0.008	54.227	55.934

Algoritmo general

Dado un estado actual \(\boldsymbol{\theta}^{(b)} = (\theta_1^{(b)},\ldots,\theta_k^{(b)})\), se genera un nuevo estado \(\boldsymbol{\theta}^{(b+1)}\) de la siguiente manera:

Muestrear \(\theta_1^{(b+1)}\sim p\left(\theta_1\mid\theta_2^{(b)},\theta_3^{(b)},\ldots,\theta_k^{(b)}\right)\).
Muestrear \(\theta_2^{(b+1)}\sim p\left(\theta_2\mid\theta_1^{(b+1)},\theta_3^{(b)},\ldots,\theta_k^{(b)}\right)\).
Muestrear \(\theta_3^{(b+1)}\sim p\left(\theta_3\mid\theta_1^{(b+1)},\theta_2^{(b+1)},\ldots,\theta_k^{(b)}\right)\).
\(\vdots\)
Muestrear \(\theta_k^{(b+1)}\sim p\left(\theta_k\mid\theta_1^{(b+1)},\theta_2^{(b+1)},\ldots,\theta_{k-1}^{(b+1)}\right)\).
Definir \(\boldsymbol{\theta}^{(b+1)} = (\theta_1^{(b+1)},\ldots,\theta_k^{(b+1)})\).
Repetir los pasos 1 a 6 hasta alcanzar la convergencia.

Cadenas de Markov

Un proceso estocástico es una colección de variables aleatorias \(\{\theta_t\in S:t\in T\}\), donde \(S\) representa el espacio de estados, que puede ser discreto, e.g., \(\{1,\ldots,k\}\), o continuo, e.g., \((-\infty,\infty)\), y \(T\) representa el conjunto de índices, que puede ser discreto, e.g., \(\{0,1,\ldots\}\), o continuo, e.g., \([0,\infty)\).

Un proceso estocástico \(\{\theta_t\in S:t\in T\}\) con \(T=\{0,1,\ldots\}\) se denomina cadena de Markov si, para todo \(A\subset S\), se cumple que
\[ \textsf{Pr}(\theta_{t+1}\in A\mid \theta_0,\ldots,\theta_t) = \textsf{Pr}(\theta_{t+1}\in A\mid\theta_t)\,. \]

Cadenas de Markov ergódicas

Una cadena de Markov ergódica posee propiedades que garantizan su convergencia hacia un comportamiento estable a largo plazo, es decir, a una distribución estacionaria.

Para que una cadena de Markov sea ergódica, debe cumplir las siguientes condiciones:

Irreducibilidad: Cualquier estado \(i\) puede alcanzarse desde cualquier otro estado \(j\) en un número finito de pasos con probabilidad positiva.
Aperiodicidad: No existe un período fijo \(d > 1\) que restrinja las visitas a cada estado a múltiplos de \(d\).
Recurrencia positiva: Todo estado \(i\) tiene un tiempo de retorno esperado finito.

Cuando estas propiedades se cumplen, la cadena de Markov converge a una única distribución estacionaria, independientemente del estado inicial.

Teorema Ergódico

Las cadenas de Markov ergódicas tienen una distribución estacionaria única \(\pi(\theta)\), la cual describe el comportamiento asintótico de la cadena tras un número suficientemente grande de iteraciones, independientemente del estado inicial.

(Teorema Ergódico.) Si una cadena de Markov \(\{\theta_t\in S:t\in T\}\) es ergódica y si \(f\) es una función real tal que \(\textsf{E}_\pi|g(\theta)|<\infty\), entonces, con probabilidad 1, cuando \(B\to\infty\), se cumple que
\[ \frac{1}{B}\sum_{b=1}^B g(\theta^{(b)}) \longrightarrow \textsf{E}_{\pi(\theta)}(g(\theta)) = \int_\Theta g(\theta)\,\pi(\theta)\,\textsf{d}\theta, \]
donde el valor esperado de \(f(\theta)\) se toma con respecto a la distribución estacionaria \(\pi(\theta)\).

Resumen

Una cadena de Markov generada mediante el muestreador de Gibbs es ergódica y tiene como distribución estacionaria la distribución posterior.

La secuencia \(\boldsymbol{\theta}^{(1)},\ldots,\boldsymbol{\theta}^{(B)}\) puede utilizarse como si fuera una muestra aleatoria de valores de \(\boldsymbol{\theta}\) extraídos de la distribución posterior \(p(\boldsymbol{\theta} \mid \boldsymbol{y})\). En particular,
\[ \frac{1}{B} \sum_{b=1}^{B} g(\boldsymbol{\theta}^{(b)}) \longrightarrow \textsf{E}[g(\boldsymbol{\theta}) \mid \boldsymbol{y}]=\int_{\Theta} g(\boldsymbol{\theta})\, p(\boldsymbol{\theta} \mid \boldsymbol{y})\,\textsf{d}\boldsymbol{\theta} \quad \text{cuando} \quad B\rightarrow \infty. \]

El punto de partida \(\boldsymbol{\theta}^{(0)}\) es arbitrario y suele elegirse a partir de la distribución previa para facilitar la convergencia.

El muestreador de Gibbs pertenece a la familia de técnicas de aproximación conocidas como cadenas de Markov de Monte Carlo (Markov Chain Monte Carlo, MCMC).

Los métodos de MCMC son técnicas de aproximación numérica, no son modelos, y no aportan información adicional más allá de la contenida en los datos \(\boldsymbol{y}\).

Diagnósticos de convergencia

El muestreador de Gibbs genera muestras que eventualmente convergen a la distribución posterior, aunque en algunos casos la convergencia puede ser lenta debido a la autocorrelación entre las muestras.

Para evaluar la calidad de la simulación, es fundamental considerar las siguientes preguntas:

¿Cuál es el punto de partida de la cadena?
¿Cuándo la cadena alcanza el equilibrio?
Una vez alcanzado el equilibrio, ¿cuánto tiempo se debe monitorear la cadena?

Dado que la convergencia no puede verificarse con certeza absoluta, solo es posible detectar cuándo la cadena aún no ha convergido.

Serie

Las series permiten evaluar la estacionariedad de la cadena.

Si se detectan problemas de estacionariedad, se recomienda aumentar el número de iteraciones para mejorar la convergencia.

# Log-verosimilitud
ll <- numeric(B)
for (i in 1:B) {
  ll[i] <- sum(dnorm(y, mean = THETA$theta[i], sd = sqrt(THETA$sigma2[i]), log = TRUE))
}

# Graficar la cadena de la log-verosimilitud
par(mfrow = c(1,1), mar = c(2.75, 3, 1.5, 0.5), mgp = c(1.7, 0.7, 0))
plot(ll[-c(1:10)], type = "p", pch = 16, col = adjustcolor(1, 0.3), cex = 0.8, 
     xlab = "Iteración", ylab = "Log-verosimilitud", 
     main = "Cadena de log-verosimilitud")

Autocorrelación

La función de autocorrelación se define como
\[ \operatorname{acf}_{t}(\theta)=\frac{\frac{1}{B-t} \sum_{b=1}^{B-t}\left(\theta^{(b)}-\bar{\theta}\right)\left(\theta^{(b+t)}-\bar{\theta}\right)}{\frac{1}{B-1} \sum_{b=1}^{B}\left(\theta^{(b)}-\bar{\theta}\right)^{2}}, \qquad\text{donde} \qquad\bar{\theta} = \frac{1}{B}\sum_{b=1}^B \theta^{(b)}. \]

Una mayor autocorrelación implica la necesidad de más muestras para alcanzar la precisión deseada.

Para reducir la autocorrelación en la cadena, se recomienda:

Adelgazar la cadena mediante muestreo sistemático.
Reparametrizar el modelo para mejorar la eficiencia de la exploración del espacio de parámetros.

# Autocorrelación
par(mfrow = c(1, 2), mar = c(3, 3, 2.5, 1), mgp = c(1.75, 0.75, 0))

# Función de autocorrelación para theta
acf(THETA$theta, main = expression("Autocorrelación de " ~ theta))

# Función de autocorrelación para sigma^2
acf(THETA$sigma2, main = expression("Autocorrelación de " ~ sigma^2))

Tamaño efectivo de muestra

Las cadenas generadas suelen presentar autocorrelación, lo que implica que aportan menos información en comparación con una secuencia de muestras independientes e idénticamente distribuidas (IID) de la distribución posterior.

El tamaño efectivo de la muestra representa la cantidad de muestras IID que contendrían la misma información que la cadena, y se define como: \[ n_{\text{eff}}(\theta) = \frac{B}{1+2\sum_{t=1}^\infty \text{acf}_t(\theta)}\,, \]
donde \(\text{acf}_t(\theta)\) es la autocorrelación en el rezago \(t\). En la práctica, la suma se trunca cuando dos términos consecutivos se vuelven negativos, es decir, cuando se cumple que
\[ \text{acf}_{t-1} + \text{acf}_t < 0. \]
Este criterio evita la inclusión de contribuciones espurias en la estimación del tamaño efectivo de la muestra.

La librería coda en R (Convergence Diagnosis and Output Analysis) se utiliza para analizar y diagnosticar la convergencia de cadenas de Markov generadas por algoritmos MCMC.

# Librería necesaria 
suppressMessages(suppressWarnings(library(coda)))

# Tamaño efectivo de la muestra
neff <- coda::effectiveSize(THETA)
round(neff, 0)

##  theta sigma2 
##  10278  10000

Error estándar de Monte Carlo

El error estándar de Monte Carlo (EMC) mide la desviación estándar del promedio muestral ajustado por el tamaño efectivo de la muestra:
\[ \textsf{EMC}(\hat\theta) = \frac{\textsf{DE}(\hat\theta)}{\sqrt{n_{\text{eff}}(\theta)}}\,, \]
donde \(\hat\theta = \frac{1}{B} \sum_{b=1}^B \theta^{(b)}\) es la media posterior, \(\textsf{DE}(\hat\theta) = \sqrt{\frac{1}{B-1} \sum_{b=1}^B (\theta^{(b)} - \hat\theta)^2}\) es la desviación estándar posterior, y \(n_{\text{eff}}\) es el tamaño efectivo de la muestra.

El EMC cuantifica la precisión de la estimación de \(\hat\theta\), disminuyendo conforme aumenta \(n_{\text{eff}\), lo que implica una mayor eficiencia en la simulación.

# Error estándar de Monte Carlo
EMC <- apply(X = THETA, MARGIN = 2, FUN = sd)/sqrt(neff)
round(EMC, 4)

##  theta sigma2 
## 0.0042 0.0736

# Coeficiente de variación de Monte Carlo (%)
CVMC <- 100*EMC/abs(colMeans(THETA))
round(CVMC, 4)

##  theta sigma2 
## 0.0077 0.0510

Prueba de Gelman-Rubin

La prueba de Gelman-Rubin evalúa la convergencia en algoritmos MCMC comparando la varianza entre cadenas con la varianza dentro de cada cadena.

Dado un conjunto de \(M\) cadenas de longitud \(B\), se define la varianza dentro de las cadenas como
\[ \textsf{W} = \frac{1}{M} \sum_{m=1}^{M} s_m^2, \quad \text{donde} \quad s_m^2 = \frac{1}{B-1} \sum_{b=1}^{B} (\theta_m^{(b)} - \bar{\theta}_m)^2. \]
y la varianza entre las cadenas como
\[ \textsf{B} = \frac{B}{M-1} \sum_{m=1}^{M} (\bar{\theta}_m - \bar{\theta})^2, \quad \text{donde} \quad \bar{\theta} = \frac{1}{M} \sum_{m=1}^{M} \bar{\theta}_m. \]

A partir de esto, el estimador de la varianza total es
\[ \hat{\sigma}^2 = \frac{B - 1}{B} \, \textsf{W} + \frac{1}{B} \, \textsf{B}, \]
y el estadístico de Gelman-Rubin se define como
\[ \hat{R} = \sqrt{\frac{\hat{\sigma}^2}{\textsf{W}}}. \]
Si \(\hat{R} \approx 1\), las cadenas han convergido; si \(\hat{R} > 1.1\), se recomienda continuar la simulación. Valores grandes indican que la exploración del espacio de parámetros es insuficiente.

gibbs_sampler <- function(y, B, mu0, t20, nu0, s20, verbose = TRUE) {
  # Cantidades derivados
  n <- length(y)
  yb <- mean(y)
  
  # Inicialización
  sig2 <- 1 / rgamma(1, shape = nu0 / 2, rate = nu0 * s20 / 2)
  THETA <- matrix(NA, nrow = B, ncol = 2)
  
  # Iteraciones del muestreador de Gibbs
  for (b in 1:B) {
    # Actualización de theta
    t2n   <- 1 / (1 / t20 + n / sig2)
    mun   <- t2n * (mu0 / t20 + n * yb / sig2)
    theta <- rnorm(1, mean = mun, sd = sqrt(t2n))
    
    # Actualización de sigma^2
    nun  <- nu0 + n 
    s2n  <- (nu0 * s20 + sum((y - theta)^2)) / nun
    sig2 <- 1 / rgamma(1, shape = nun / 2, rate = nun * s2n / 2)
    
    # Almacenar muestras
    THETA[b, ] <- c(theta, sig2)
    
    # Mostrar progreso cada 10% de las iteraciones si verbose = TRUE
    if (verbose && (b %% floor(B / 10) == 0)) 
      cat(sprintf("%.0f%% completado...\n", 100 * b / B))
  }
  
  # Convertir THETA en data frame y asignar nombres a las columnas
  THETA <- as.data.frame(THETA)
  colnames(THETA) <- c("theta", "sigma2")
  
  return(THETA)
}

# Cadenas
set.seed(123)
cadena_1 <- gibbs_sampler(y, B, mu0, t20, nu0, s20, verbose = F)
cadena_2 <- gibbs_sampler(y, B, mu0, t20, nu0, s20, verbose = F)
cadena_3 <- gibbs_sampler(y, B, mu0, t20, nu0, s20, verbose = F)

# Crear un objeto mcmc.list con múltiples cadenas
chains_mcmc <- mcmc.list(
     mcmc(data = cadena_1, start = 10, end = B, thin = 1), 
     mcmc(data = cadena_2, start = 10, end = B, thin = 1), 
     mcmc(data = cadena_3, start = 10, end = B, thin = 1)
)

# Calcular R-hat con coda
gelman.diag(chains_mcmc)$psrf

##        Point est. Upper C.I.
## theta   1.0001892   1.000867
## sigma2  0.9999682   1.000075

Más diagnósticos

Los siguientes diagnósticos también se encuentran disponibles en la librería coda de R:

heidel.diag() : Prueba estacionariedad y recomienda descartar un número óptimo de iteraciones iniciales para asegurar convergencia.
geweke.diag() : Compara medias al inicio y al final de la cadena para detectar diferencias significativas que indiquen falta de convergencia.
raftery.diag(): Estima el número mínimo de iteraciones necesarias para garantizar intervalos de credibilidad con una precisión y cobertura deseada.

Ejercicios

Demuestre que la distribución t se puede expresar como una mezcla de distribuciones Normales ponderadas por una distribución Gamma-Inversa. Es decir, demostrar que la distribución muestral \(y_{i}\mid\theta,\sigma^2 \stackrel{\text {iid}}{\sim} \textsf{t}_\kappa(\theta,\sigma^2)\), para \(i=1,\ldots,n\), es equivalente a la distribución jerárquica dada por \[ y_{i}\mid\theta,V_i \stackrel{\text {ind}}{\sim} \textsf{N}(\theta,V_i)\,,\qquad V_{i}\mid\sigma^2\stackrel{\text {iid}}{\sim} \textsf{GI}\left( \frac{\kappa}{2}, \frac{\kappa\,\sigma^2}{2}\right)\,. \]
Considere los datos sobre el número de hijos de hombres en sus 30s, con y sin título universitario, reportados en menchild30bach.dat y menchild30nobach.dat. Asuma modelos Poisson para ambos grupos, pero ahora parametrizamos \(\theta_A\) y \(\theta_B\) como \(\theta_A = \theta\) y \(\theta_B = \theta \times \phi\), donde \(\phi\) representa la razón relativa \(\theta_B / \theta_A\). Suponemos distribuciones previas \(\theta \sim \textsf{Gamma}(a_\theta, b_\theta)\) y \(\phi \sim \textsf{Gamma}(a_\phi, b_\phi)\).
1. ¿Son \(\theta_A\) y \(\theta_B\) independientes o dependientes bajo esta distribución previa? ¿En qué situaciones se justifica el uso de una distribución previa de este tipo?
2. Obtenga la distribución condicional completa de \(\theta\).
3. Obtenga la distribución condicional completa de \(\phi\).
4. Fije \(a_\theta = 2\) y \(b_\theta = 1\). Considere \(a_\phi = b_\phi \in \{8, 16, 32, 64, 128\}\). Para cada uno de estos cinco valores, implemente un muestreador de Gibbs con al menos 10,000 iteraciones. Realice un análisis exhaustivo de convergencia.
5. Estime \(\textsf{E}(\theta_B - \theta_A \mid \boldsymbol{y}_A, \boldsymbol{y}_B)\) en cada caso. Describa los efectos de la distribución previa de \(\phi\) sobre los resultados
El archivo glucose.dat contiene las concentraciones de glucosa en plasma de 532 mujeres en un estudio sobre diabetes.
1. Realice un histograma de los datos. Describa cómo la distribución empírica difiere de una distribución Normal en términos de asimetría, curtosis u otras características relevantes.
2. Considere el siguiente modelo de mezcla para estos datos. Cada participante del estudio tiene asociada una variable latente (no observada) de pertenencia a un grupo \(\xi_i\), que sigue una distribución Bernoulli con probabilidad \(\omega\), de forma que \(\xi_i = 1\) con probabilidad \(\omega\) y \(\xi_i = 2\) con probabilidad \(1 - \omega\). Si \(\xi_i = 1\), entonces \(y_i \mid \theta_1, \sigma_1^2 \sim \textsf{N}(\theta_1, \sigma_1^2)\), y si \(\xi_i = 2\), entonces \(y_i \mid \theta_2, \sigma_2^2 \sim \textsf{N}(\theta_2, \sigma_2^2)\). Las distribuciones previas están dadas por \(\omega \sim \textsf{Beta}(\alpha_0, \beta_0)\), \(\theta_j \sim \textsf{N}(\mu_0, \tau_0^2)\), \(\sigma_j^2 \sim \textsf{GI}(\nu_0/2, \nu_0\,\sigma_0^2/2)\), para \(j = 1,2\). Obtenga las distribuciones condicionales completas de \(\xi_1, \dots, \xi_n\), \(\omega\), \(\theta_1\), \(\theta_2\), \(\sigma_1^2\) y \(\sigma_2^2\).
3. Fijando \(\alpha_0 = \beta_0 = 1\), \(\mu_0 = 120\), \(\tau_0^2 = 200\), \(\sigma_0^2 = 1000\) y \(\nu_0 = 10\), implemente un muestreador de Gibbs con al menos 10,000 iteraciones. Realice un análisis exhaustivo de convergencia.
4. Para cada iteración \(b\) del muestreador de Gibbs, denere un valor \(\xi \sim \textsf{Ber}(\omega^{(b)})\) y luego muestree \(y^{*\,(b)} \sim \textsf{N}(\theta_{\xi}^{(b)}, \sigma_{\xi}^{2\,(b)})\), para \(b = 1,\ldots,B\). Por medio de un histograma, represente la distribución de \(y^{*\,(1)}, \dots, y^{*\,(B)}\) y compárela con la distribución obtenida en la parte a. Discuta la adecuación de este modelo de mezcla de dos componentes para los datos de glucosa.
Un estudio de panel siguió a 25 parejas casadas durante cinco años para analizar la relación entre las tasas de divorcio y diversas características de las parejas. Los investigadores buscan modelar la probabilidad de divorcio en función de estas características. Los datos están disponibles en el archivo divorce.dat. Se utilizará regresión probit, donde una variable binaria \(y_i\) se modela mediante la siguiente formulación usando variables latentes:
\[ z_i = \beta x_i + \epsilon_i, \quad y_i = I(z_i > \delta), \quad \text{para } i = 1,\ldots,n, \] donde \(\beta\) y \(\delta\) son parámetros desconocidos, \(\epsilon_i \overset{\text{iid}}{\sim} \textsf{N}(0,1)\) para \(i=1,\ldots,n\), y \(I(z_i > \delta)\) es una función indicadora que toma el valor 1 si \(z_i > \delta\) y 0 en caso contrario.
1. Determine la distribución condicional completa de \(\beta\), asumiendo que \(\beta \sim \textsf{N}(0, \tau^2_0)\).
2. Determine que la distribución condicional completa de \(\delta\), asumiendo que \(\delta \sim \textsf{N}(0, \kappa^2_0)\).
3. Determine la distribución condicional completa de \(z_i\).
4. Fijando \(\tau^2_0 = \kappa^2_0 = 16\), implemente un muestreador de Gibbs para aproximar la distribución posterior conjunta de \(z_1,\ldots z_n\), \(\beta\) y \(\delta\). Ejecute el muestreador de Gibbs el tiempo suficiente para que los tamaños efectivos de muestra de todos los parámetros desconocidos, incluidos los \(z_i\), sean mayores a 1,000. Realice un análisis exhaustivo de convergencia.
5. Obtenga un intervalo de credibilidad posterior del 95% para \(\beta\), así como \(\textsf{Pr}(\beta > 0 \mid \boldsymbol{x}, \boldsymbol{y})\). Discuta los resultados obtenidos.
Sea \(y_t = \rho y_{t-1} + \epsilon_t\), donde \(\epsilon_t\) es un proceso i.i.d. con \(\epsilon_t \sim \textsf{N}(0, \sigma^2)\). Este es un modelo ampliamente utilizado en el análisis de series temporales conocido como el modelo autorregresivo de orden uno (AR(1)).
1. Escriba la verosimilitud condicional dada \(y_1\), es decir, \(p(y_2, \dots, y_n \mid y_1, \rho, \sigma^2)\).
2. Suponga una distribución a priori de la forma \(p(\rho, \sigma^2) \propto 1/\sigma^2\):
  - Determine la distribución posterior conjunta \(p(\rho, \sigma^2 \mid y_1, \dots, y_n)\) basada en la verosimilitud condicional.
  - Obtenga \(p(\rho \mid \sigma^2, y_1, \dots, y_n)\) y \(p(\sigma^2 \mid y_1, \dots, y_n)\) utilizando la verosimilitud condicional.
  - Simule dos conjuntos de datos con \(n = 500\) observaciones cada uno: uno con \(\rho = 0.95\), \(\sigma^2 = 4\) y otro con \(\rho = 0.3\), \(\sigma^2 = 4\). Ajuste el modelo a ambos conjuntos de datos y resuma los resultados posteriores en cada caso.
La base de datos personas.csv, disponible en la página web del curso, es una muestra del módulo de Personas de la encuesta Medición de Pobreza Monetaria y Desigualdad 2021, realizada por el DANE en Colombia (enlace oficial). La muestra incluye a la población civil no institucional residente en el país, con datos recolectados mediante informantes directos (mayores de 18 años o menores que trabajen) o informantes idóneos del hogar. Se considera el ingreso total (ingtot), que corresponde a la suma de todas las fuentes de ingresos, tanto observadas como imputadas, específicamente para los habitantes de Bogotá.
1. Para ajustar el modelo \[ y_i \mid \theta, \sigma^2 \overset{\text{iid}}{\sim} \textsf{t}_\kappa(\theta, \sigma^2), \quad \text{para } i = 1, \dots, n, \] se utilizan las distribuciones previas \[ \theta \sim \textsf{N}(\mu_0, \gamma_0^2), \quad \sigma^2 \sim \textsf{Gamma}\left(\frac{\alpha_0}{2}, \frac{\beta_0}{2}\right), \quad \kappa \sim \textsf{U}\{1, 2, \dots, \nu_0\}, \] donde \(\mu_0\), \(\gamma_0^2\), \(\alpha_0\), \(\beta_0\) y \(\nu_0\) son hiperparámetros del modelo. La distribución muestral \(y_i \mid \theta, \sigma^2 \overset{\text{iid}}{\sim} \textsf{t}_\kappa(\theta, \sigma^2)\) se puede reformular de manera jerárquica como
  \[ y_i \mid \theta, \zeta_i^2 \overset{\text{ind}}{\sim} \textsf{N}(\theta, \zeta_i^2), \quad \zeta_i^2 \mid \sigma^2 \overset{\text{iid}}{\sim} \textsf{GI}\left(\frac{\kappa}{2}, \frac{\kappa \sigma^2}{2} \right), \]
  donde las variables auxiliares \(\zeta_i^2\), aunque desconocidas, se introducen para facilitar la implementación del muestreador de Gibbs. Su incorporación permite que las distribuciones condicionales completas de los parámetros desconocidos, incluidas las auxiliares, tengan formas probabilísticas estándar, lo que simplifica significativamente el proceso de muestreo.
2. Realizar inferencia sobre la media, la desviación estándar y el coeficiente de variación de los ingresos en escala logarítmica. Dado que, si \(X \mid \kappa, \theta, \sigma^2 \sim \textsf{t}_\kappa(\theta, \sigma^2)\), entonces
  \[ \textsf{E}(X) = \theta, \quad \text{para } \kappa > 1, \qquad \textsf{Var}(X) = \frac{\kappa}{\kappa - 2} \sigma^2, \quad \text{para } \kappa > 2, \]
  se puede calcular la desviación estándar como \(\textsf{DE}(X) = \sqrt{\textsf{Var}(X)}\) y el coeficiente de variación como \(\textsf{CV}(X) = \textsf{DE}(X) / \textsf{E}(X)\).
3. En el modelo jerárquico, las variables auxiliares \(\zeta_i^2\) juegan un papel fundamental en la detección de valores atípicos, ya que representan la varianza específica de cada observación \(y_i\). Estas variables permiten que el modelo ajuste dinámicamente la dispersión en torno a la media \(\theta\), asignando una mayor varianza a las observaciones que se desvían significativamente del patrón general. Este mecanismo mejora la identificación de outliers sin distorsionar las estimaciones globales del modelo. Aplique los siguientes métodos para detectar valores atípicos en los ingresos (en escala logarítmica):
  - Criterio basado en \(\zeta_i^2\): Identifique las observaciones asociadas con valores de \(\zeta_i^2\) considerablemente altos en comparación con el resto. Como referencia, se pueden considerar outliers aquellas observaciones cuyo valor medio posterior de \(\zeta_i^2\) supere un umbral, como el percentil 95 de la distribución estimada de \(\zeta_i^2\). Estas observaciones son indicativas de outliers, ya que el modelo les asigna una varianza elevada para reflejar su desviación respecto a la media \(\theta\).
  - Criterio basado en residuos estandarizados: Calcule los residuos estandarizados utilizando la fórmula: \[ r_i = \frac{y_i - \theta}{\sqrt{\zeta_i^2}}, \] donde valores absolutos grandes de \(r_i\), típicamente \(|r_i| > 3\), sugieren la presencia de outliers. Este enfoque combina la desviación de \(y_i\) respecto a \(\theta\) y la varianza específica \(\zeta_i^2\), proporcionando una medida confiable para identificar valores extremos.
  - Ambos criterios pueden ser complementados mediante visualizaciones, como gráficos de los valores estimados de \(\zeta_i^2\) o de los residuos estandarizados \(r_i\). Estas herramientas permiten una exploración más intuitiva y detallada de los posibles outliers en el conjunto de datos.

Referencias

Hoff, P. D. (2009). A First Course in Bayesian Statistical Methods. Springer New York.

Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2013). Bayesian Data Analysis (3rd ed.). Chapman & Hall/CRC.

Muestreador de Gibbs

Juan Sosa PhD

Email jcsosam@unal.edu.co

GitHub https://github.com/jstats1702

Introducción

Muetreador de Gibbs

Distribuciones condicionales completas

Ejemplo: Puntajes de Matemáticas

Tratamiento de datos

Muestreador de Gibbs

Cadenas

Distribución posterior

Distriución predictiva

Chequeo del modelo

Inferencia: Bogotá

Algoritmo general

Cadenas de Markov

Cadenas de Markov ergódicas

Teorema Ergódico

Resumen

Diagnósticos de convergencia

Serie

Autocorrelación

Tamaño efectivo de muestra

Error estándar de Monte Carlo

Prueba de Gelman-Rubin

Más diagnósticos

Ejercicios

Referencias