Modelo Normal Multivariado

Introducción

El modelo Normal multivariado para variables multivariadas \(\boldsymbol{y}_1, \ldots, \boldsymbol{y}_n\), donde \(\boldsymbol{y}_i = (y_{i,1}, \ldots, y_{i,p}) \in \mathbb{R}^p\), con \(i = 1, \ldots, n\), se define como: \[ \begin{aligned} \boldsymbol{y}_i \mid \boldsymbol{\theta}, \mathbf{\Sigma} &\stackrel{\text{iid}}{\sim} \textsf{N}(\boldsymbol{\theta}, \mathbf{\Sigma}) \\ (\boldsymbol{\theta}, \mathbf{\Sigma}) &\sim p(\boldsymbol{\theta}, \mathbf{\Sigma}) \end{aligned} \] donde \(\boldsymbol{\theta}\) es el vector de medias y \(\mathbf{\Sigma}\) es la matriz de covarianzas.

El modelo involucra \(k = p + \frac{p(p+1)}{2}\) parámetros desconocidos por estimar.

El vector aleatorio \(\boldsymbol{y} = (y_1, \ldots, y_p)\) sigue una distribución Normal multivariada si su función de densidad de probabilidad está dada por: \[ p(\boldsymbol{y} \mid \boldsymbol{\theta}, \mathbf{\Sigma}) = (2\pi)^{-p/2} |\mathbf{\Sigma}|^{-1/2} \exp\left\{ -\frac{1}{2} (\boldsymbol{y} - \boldsymbol{\theta})^{\textsf{T}} \mathbf{\Sigma}^{-1} (\boldsymbol{y} - \boldsymbol{\theta}) \right\}, \] donde: \[ \boldsymbol{\theta} = (\theta_1, \ldots, \theta_p), \qquad \mathbf{\Sigma} = \begin{bmatrix} \sigma_1^2 & \sigma_{1,2} & \cdots & \sigma_{1,p} \\ \sigma_{2,1} & \sigma_2^2 & \cdots & \sigma_{2,p} \\ \vdots & \vdots & \ddots & \vdots \\ \sigma_{p,1} & \sigma_{p,2} & \cdots & \sigma_p^2 \\ \end{bmatrix}. \]

La matriz \(\mathbf{\Sigma}\) es simétrica, \(\mathbf{\Sigma}^{\textsf{T}} = \mathbf{\Sigma}\), y definida positiva, \(\boldsymbol{x}^{\textsf{T}} \mathbf{\Sigma} \boldsymbol{x} > 0\), para todo \(\boldsymbol{x} \in \mathbb{R}^p\).

(Ejercicio.) El núcleo de la distribución Normal multivariada está dado por: \[ p(\boldsymbol{y} \mid \boldsymbol{\theta}, \mathbf{\Sigma}) \propto \exp\left\{ -\frac{1}{2} \left[ \boldsymbol{y}^{\textsf{T}} \mathbf{\Sigma}^{-1} \boldsymbol{y} - 2 \boldsymbol{y}^{\textsf{T}} \mathbf{\Sigma}^{-1} \boldsymbol{\theta} \right] \right\}. \]

Estadístico suficiente

(Ejercicio.) Si \(\boldsymbol{y}_i \mid \boldsymbol{\theta}, \mathbf{\Sigma} \stackrel{\text{iid}}{\sim} \textsf{N}(\boldsymbol{\theta}, \mathbf{\Sigma})\), con \(i = 1, \ldots, n\), entonces la distribución muestral conjunta de las observaciones está dada por \[ p\left(\mathbf{Y} \mid \boldsymbol{\theta}, \mathbf{\Sigma} \right) = \left(2 \pi\right)^{-np / 2} |\mathbf{\Sigma}|^{-n / 2} \exp\left\{ -\frac{1}{2} \sum_{i=1}^n (\boldsymbol{y}_i - \boldsymbol{\theta})^{\textsf{T}} \mathbf{\Sigma}^{-1} (\boldsymbol{y}_i - \boldsymbol{\theta}) \right\}, \] donde \(\mathbf{Y} = [\boldsymbol{y}_1^{\textsf{T}}, \ldots, \boldsymbol{y}_n^{\textsf{T}}]^{\textsf{T}}\).

(Ejercicio.) La expresión para la suma cuadrática en el exponente es equivalente a \[ \sum_{i=1}^n (\boldsymbol{y}_i - \boldsymbol{\theta})^{\textsf{T}} \mathbf{\Sigma}^{-1} (\boldsymbol{y}_i - \boldsymbol{\theta}) = \textsf{tr}\left( \mathbf{\Sigma}^{-1} \sum_{i=1}^n \boldsymbol{y}_i \boldsymbol{y}_i^{\textsf{T}} \right) - 2 \left(\sum_{i=1}^n \boldsymbol{y}_i\right)^{\textsf{T}} \mathbf{\Sigma}^{-1} \boldsymbol{\theta} + n \boldsymbol{\theta}^{\textsf{T}} \mathbf{\Sigma}^{-1} \boldsymbol{\theta}\,, \] lo cual sugiere que \[ \left( \sum_{i=1}^n \boldsymbol{y}_i, \sum_{i=1}^n \boldsymbol{y}_i \boldsymbol{y}_i^{\textsf{T}} \right) \] es un estadístico suficiente para \((\boldsymbol{\theta}, \mathbf{\Sigma})\).

(Ejercicio.) Además, la media muestral \(\bar{\boldsymbol{y}}\) y la matriz de covarianza muestral \(\mathbf{S}\) también forman un estadístico suficiente para \((\boldsymbol{\theta}, \mathbf{\Sigma})\), dado que: \[ \sum_{i=1}^n (\boldsymbol{y}_i - \boldsymbol{\theta})^{\textsf{T}} \mathbf{\Sigma}^{-1} (\boldsymbol{y}_i - \boldsymbol{\theta}) = \textsf{tr}\left( \mathbf{\Sigma}^{-1} \left[ (n-1) \mathbf{S} + n (\bar{\boldsymbol{y}} - \boldsymbol{\theta})(\bar{\boldsymbol{y}} - \boldsymbol{\theta})^{\textsf{T}} \right] \right), \] donde: \[ \bar{\boldsymbol{y}} = \frac{1}{n} \sum_{i=1}^n \boldsymbol{y}_i \quad \text{y} \quad \mathbf{S} = \frac{1}{n-1} \sum_{i=1}^n (\boldsymbol{y}_i - \bar{\boldsymbol{y}})(\boldsymbol{y}_i - \bar{\boldsymbol{y}})^{\textsf{T}}. \]

Modelo Normal multivariado

Para modelar una colección de observaciones normales multivariadas bajo un enfoque Bayesiano semiconjugado, se especifica el estado de información previo sobre \(\boldsymbol{\theta}\) de manera independiente de \(\mathbf{\Sigma}\), de modo que \(p(\boldsymbol{\theta},\mathbf{\Sigma})=p(\boldsymbol{\theta})\,p(\mathbf{\Sigma})\), donde \[ \boldsymbol{\theta} \sim \textsf{N}(\boldsymbol{\mu}_0, \mathbf{\Lambda}_0)\qquad\text{y}\qquad\mathbf{\Sigma} \sim \textsf{WI}(\nu_0, \mathbf{S}_0^{-1})\,, \] siendo \(\boldsymbol{\mu}_0\), \(\mathbf{\Lambda}_0\), \(\nu_0\) y \(\mathbf{S}_0\) los hiperparámetros del modelo.

Distribución Wishart

La matriz aleatoria \(\mathbf{W} > 0\) de dimensión \(p \times p\) sigue una distribución Wishart con parámetros \(\nu > p - 1\) (grados de libertad) y \(\mathbf{S} > 0\) (matriz de escala de \(p \times p\)), denotada como \(\mathbf{W} \sim \textsf{W}(\nu, \mathbf{S})\), si su función de densidad de probabilidad está dada por: \[ p(\mathbf{W} \mid \nu, \mathbf{S}) \propto |\mathbf{W}|^{(\nu - p - 1) / 2} \exp\left\{ -\frac{1}{2} \textsf{tr}(\mathbf{S}^{-1} \mathbf{W}) \right\}. \]

(Ejercicio.) Si \(\mathbf{W} \sim \textsf{W}(\nu, \mathbf{S})\), entonces la matriz de esperanza está dada por: \[ \textsf{E}(\mathbf{W}) = \nu\,\mathbf{S}. \]

Distribución Wishart Inversa

La matriz aleatoria \(\mathbf{W} > 0\) de dimensión \(p \times p\) sigue una distribución Wishart Inversa con parámetros \(\nu > p+1\) (grados de libertad) y \(\mathbf{S} > 0\) (matriz de escala de \(p \times p\)), denotada como \(\mathbf{W} \sim \textsf{WI}(\nu, \mathbf{S})\), si su función de densidad de probabilidad está dada por: \[ p(\mathbf{W} \mid \nu, \mathbf{S}^{-1}) \propto |\mathbf{W}|^{-(\nu + p + 1)/2} \exp\left\{ -\frac{1}{2} \textsf{tr}(\mathbf{S} \mathbf{W}^{-1}) \right\}. \]

(Ejercicio.) Si \(\mathbf{W} \sim \textsf{WI}(\nu, \mathbf{S}^{-1})\), entonces su valor esperado es: \[ \textsf{E}(\mathbf{W}) = \frac{1}{\nu - p - 1} \mathbf{S}, \qquad \text{para } \nu > p + 1. \]

(Ejercicio.) Si \(\mathbf{W} \sim \textsf{W}(\nu, \mathbf{S})\), entonces \(\mathbf{W}^{-1} \sim \textsf{WI}(\nu, \mathbf{S}^{-1})\).

Entrenamiento

El modelo se puede entrenar por medio de un muestreador de Gibbs (Gibbs sampler) que permita generar muestras correlacionadas de la distribución posterior por medio de las distribuciones condicionales completas de los parámetros.

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

1. Muestrear \(\boldsymbol{\theta}^{(b+1)} \sim p(\boldsymbol{\theta} \mid \mathbf{\Sigma}^{(b)}, \mathbf{Y})\).
   
2. Muestrear \(\mathbf{\Sigma}^{(b+1)} \sim p(\mathbf{\Sigma} \mid \boldsymbol{\theta}^{(b+1)}, \mathbf{Y})\).
3. Actualizar \(\Theta^{(b+1)} = (\boldsymbol{\theta}^{(b+1)}, \mathbf{\Sigma}^{(b+1)})\).
   
4. Repetir los pasos 1 a 3 hasta alcanzar convergencia.

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

A continuación se detallan las distribuciones condicionales completas para un modelo normal multivariado:

• (Ejercicio.) La distribución condicional completa de \(\boldsymbol{\theta}\) es \(\boldsymbol{\theta} \mid \mathbf{\Sigma}, \mathbf{Y} \sim \textsf{N}(\boldsymbol{\mu}_n, \mathbf{\Lambda}_n)\), donde:
  \[ \boldsymbol{\mu}_n = \left( \mathbf{\Lambda}_0^{-1} + n \mathbf{\Sigma}^{-1} \right)^{-1} \left( \mathbf{\Lambda}_0^{-1} \boldsymbol{\mu}_0 + n \mathbf{\Sigma}^{-1} \bar{\boldsymbol{y}} \right) \qquad\text{y}\qquad \mathbf{\Lambda}_n = \left( \mathbf{\Lambda}_0^{-1} + n \mathbf{\Sigma}^{-1} \right)^{-1}. \]

• (Ejercicio.) La distribución condicional completa de \(\mathbf{\Sigma}\) es \(\mathbf{\Sigma} \mid \boldsymbol{\theta}, \mathbf{Y} \sim \textsf{WI}(\nu_n, \mathbf{S}_n^{-1})\), donde:
  \[ \nu_n = \nu_0 + n, \qquad \mathbf{S}_n = \mathbf{S}_0 + \sum_{i=1}^n (\boldsymbol{y}_i - \boldsymbol{\theta})(\boldsymbol{y}_i - \boldsymbol{\theta})^{\textsf{T}}. \]

(Ejercicio.) El término \(\mathbf{S}_n\) puede reescribirse en términos de la matriz de covarianza muestral \(\mathbf{S}\) y la media muestral \(\bar{\boldsymbol{y}}\) como:
\[ \mathbf{S}_n = \mathbf{S}_0 + (n-1) \mathbf{S} + n (\bar{\boldsymbol{y}} - \boldsymbol{\theta})(\bar{\boldsymbol{y}} - \boldsymbol{\theta})^{\textsf{T}}. \]

Estas distribuciones condicionales completas son la base para actualizar secuencialmente los parámetros \(\boldsymbol{\theta}\) y \(\mathbf{\Sigma}\) en cada iteración del algoritmo.

Ejemplo: Comprensión de lectura

Se realiza un estudio con una muestra de 22 niños a quienes se aplican pruebas de comprensión lectora antes y después de recibir un método de instrucción específico.

Para cada estudiante \(i\), se registran dos variables: \(y_{i,1}\) y \(y_{i,2}\), que representan los puntajes obtenidos antes y después de la instrucción, respectivamente.

El objetivo principal del análisis es evaluar la efectividad del método de enseñanza y examinar la consistencia de la prueba de comprensión lectora, es decir, determinar si los puntajes reflejan cambios significativos y coherentes asociados al método empleado.

Referencia: Hoff, P. D. (2009). A first course in Bayesian statistical methods (Vol. 580). New York: Springer.

Tratamiento de datos

# Datos: puntajes de comprensión de lectura antes y después de la instrucción
Y <- matrix(
  c(
    59, 43, 34, 32, 42, 38, 55, 67, 64, 45, 49, 72, 34, 70, 34, 50, 41, 52, 60, 34, 28, 35,
    77, 39, 46, 26, 38, 43, 68, 86, 77, 60, 50, 59, 38, 48, 55, 58, 54, 60, 75, 47, 48, 33
  ),
  nrow = 22, ncol = 2,
  dimnames = list(NULL, c("pre_test", "post_test"))
)

# Dimensiones
(n <- nrow(Y))

## [1] 22

(p <- ncol(Y))

## [1] 2

# Inspeccionar datos
summary(Y)

##     pre_test       post_test    
##  Min.   :28.00   Min.   :26.00  
##  1st Qu.:34.25   1st Qu.:43.75  
##  Median :44.00   Median :52.00  
##  Mean   :47.18   Mean   :53.86  
##  3rd Qu.:58.00   3rd Qu.:60.00  
##  Max.   :72.00   Max.   :86.00

# Estadísticos suficientes
yb <- c(colMeans(Y))
round(yb, 1)

##  pre_test post_test 
##      47.2      53.9

SS <- cov(Y)
round(SS, 1) 

##           pre_test post_test
## pre_test     182.2     148.4
## post_test    148.4     243.6

Elicitación de lo hiperparámetros

El examen fue diseñado para otorgar puntajes con un promedio de 50 sobre 100, lo que define el vector de medias previas como \(\boldsymbol{\mu}_0 = (50, 50)\).

Se utiliza una varianza previa que asegura que \(\textsf{Pr}(0 < \theta_j < 100) \approx 0.99\). Esto se logra fijando \(\sigma^2_{0,1} = \sigma^2_{0,2} = \left(\frac{50}{3}\right)^2 \approx 278\).

Adicionalmente, se establece una correlación previa de \(\rho_0 = 0.5\), lo que implica que la covarianza previa es \(\sigma_{0,12} = (0.5) \left(\frac{50}{3}\right)^2 \approx 139\).

\(\nu_0 = 4\) se elige porque es el mínimo valor entero requerido para garantizar que \(\textsf{E}(\mathbf{\Sigma})\) existe bajo una distribución Wishart Inversa.

\(\mathbf{S}_0 = \mathbf{\Lambda}_0\) asegura consistencia con las suposiciones previas sobre \(\mathbf{\Sigma}\).

# Previa
mu0 <- c(50,50)
L0  <- matrix(data = c(278,139,139,278), nrow = 2, ncol = 2)
nu0 <- 4
S0  <- matrix(data = c(278,139,139,278), nrow = 2, ncol = 2)

Ajuste del Modelo Normal Multivariado

# Inicialización
Sigma <- solve(rWishart::rWishart(n = 1, df = nu0, Sigma = S0)[,,1])

# Número de muestras
B <- 10000

# Almacenamiento
THETA <- NULL
SIGMA <- NULL
YS    <- NULL
LL    <- numeric(B)

# Cálculo de cantidades constantes previas l procedimiento iterativo
iL0 <- solve(L0)
Lm0 <- iL0 %*% mu0
nun <- nu0 + n
SSn <- S0 + (n - 1) * SS

# Muestreador de Gibbs
set.seed(1234)
for (b in 1:B) {
  # Actualización de theta
  iSigma <- solve(Sigma)
  Ln     <- solve(iL0 + n * iSigma)
  theta  <- c(mvtnorm::rmvnorm(n = 1, mean = Ln %*% (Lm0 + n * (iSigma %*% yb)), sigma = Ln))

  # Actualización de Sigma
  Sigma <- MCMCpack::riwish(v = nun, S = SSn + n * ((yb - theta) %*% t(yb - theta)))

  # Muestra de la predictiva posterior
  YS <- rbind(YS, mvtnorm::rmvnorm(n = 1, mean = theta, sigma = Sigma))

  # Cálculo de log-verosimilitud
  LL[b] <- sum(apply(X = Y, MARGIN = 1, FUN = function(y) mvtnorm::dmvnorm(y, mean = theta, sigma = Sigma, log = T)))

  # Almacenamiento de resultados
  THETA <- rbind(THETA, theta)
  SIGMA <- rbind(SIGMA, c(Sigma))

  # Mostrar progreso cada 10%
  if (b %% floor(B / 10) == 0) 
    cat(100 * round(b / B, 1), "% completado...\n", sep = "")
}

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

# Etiquetado de columnas para claridad
colnames(THETA) <- c("theta1", "theta2")
colnames(SIGMA) <- c("sigma1^2", "sigma21", "sigma12", "sigma2^2")

Convergencia

Se evalúa la convergencia del modelo por medio de los tamaños efectivos de muestra, los errores estándar y los coeficientes de variación:

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

##   theta1   theta2 sigma1^2  sigma21  sigma12 sigma2^2 
##    10000    10000     9059     9013     9013     8951

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

##   theta1   theta2 sigma1^2  sigma21  sigma12 sigma2^2 
##   0.0284   0.0322   0.6120   0.6007   0.6007   0.8182

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

##   theta1   theta2 sigma1^2  sigma21  sigma12 sigma2^2 
##   0.0603   0.0601   0.3295   0.4076   0.4076   0.3340

Inferencia

Distribución posterior de \(\boldsymbol{\theta}\) y distribución predictiva posterior:

Inferencia sobre \(\theta_2 - \theta_1\):

¿Cuál es la probabilidad posterior de que la calificación promedio del segundo examen sea mayor que la del primero?

# Diferencias theta2 - theta1
round(mean(THETA[,2] - THETA[,1] > 0), 4)

## [1] 0.9932

# Cuantiles de la diferencia posterior
round(quantile(THETA[,2] - THETA[,1], probs = c(0.025, 0.5, 0.975)), 4)

##    2.5%     50%   97.5% 
##  1.7158  6.5206 11.3320

Inferencia sobre \(\tilde{y}_2 - \tilde{y}_1\):

¿Cuál es la probabilidad posterior de que un niño seleccionado al azar obtenga una puntuación más alta en el segundo examen que en el primero?

# Probabilidad de mejora (y2 > y1)
round(mean(YS[,2] - YS[,1] > 0), 4)

## [1] 0.7121

# Cuantiles de la diferencia predictiva
round(quantile(YS[,2] - YS[,1], probs = c(0.025, 0.5, 0.975)), 4)

##     2.5%      50%    97.5% 
## -16.5155   6.4438  29.7073

Inferencia sobre \(\rho=\frac{\sigma_{1,2}}{\sigma_{1}\,\sigma_{2}}\):

¿Las pruebas son consistentes? ¿Cuál es la probabilidad posterior de que la correlación entre las calificaciones sea superior a 0.6?

# Muestras de rho
RHO <- SIGMA[,2]/sqrt(SIGMA[,1]*SIGMA[,4])

# Probabilidad posterior de que rho > 0.5
round(mean(RHO > 0.5), 4)

## [1] 0.9384

# Cuantiles de rho (2.5%, 50%, 97.5%)
round(quantile(RHO, prob = c(0.025,0.5,0.975)), 4)

##   2.5%    50%  97.5% 
## 0.4238 0.6969 0.8542

Bondad de ajuste

Se evalúa la bondad de ajuste del modelo tomando como estadísticos de prueba la diferencia de medias y el coeficiente de correlación:

# Matriz para almacenar los estadísticos
TS <- matrix(NA, B, 2)

set.seed(123)
for (b in 1:B) {
  # Generación de datos replicados
  Yrep <- mvtnorm::rmvnorm(
       n = n, 
       mean = THETA[b, ], 
       sigma = matrix(SIGMA[b, ], nrow = p, ncol = p)
  )
  
  # Cálculo de los estadísticos de interés
  TS[b, ] <- c(mean(Yrep[, 2]) - mean(Yrep[, 1]), cor(Yrep[, 1], Yrep[, 2]))
}

# valores ppp
(ppp <- c(
     mean(TS[, 1] < (mean(Y[, 2] - Y[, 1]))), 
     mean(TS[, 2] < cor(Y[, 1], Y[, 2]))
))

## [1] 0.5171 0.4976

Ejercicios

• La distribución previa de Jeffreys para el modelo Normal multivariado es:
  \[ p_J(\boldsymbol{\theta}, \mathbf{\Sigma}) \propto |\mathbf{\Sigma}|^{-(p+2)/2}. \]

  1. Explique por qué la función \(p_J\) no es una densidad de probabilidad válida para \((\boldsymbol{\theta}, \mathbf{\Sigma})\).
  2. Determine la distribución posterior de \(\boldsymbol{\theta}\) y \(\mathbf{\Sigma}\).
  3. Determine la distribución condicional completa de \(\boldsymbol{\theta}\).
  4. Determine la distribución condicional completa de \(\mathbf{\Sigma}\).

• Sea \(\mathbf{\Omega} = \mathbf{\Sigma}^{-1}\) la matriz de precisión del modelo Normal multivariado.

  1. Demuestre que una distribución previa de información unitaria para \((\boldsymbol{\theta}, \mathbf{\Omega})\) está dada por \(\boldsymbol{\theta} \mid \mathbf{\Omega} \sim \textsf{N}(\bar{\boldsymbol{y}}, \mathbf{\Omega}^{-1})\) y \(\mathbf{\Omega} \sim \textsf{W}(p+1, \mathbf{S}^{-1})\), donde \(\mathbf{S} = \frac{1}{n} \sum_{i=1}^{n} (\boldsymbol{y}_i - \bar{\boldsymbol{y}})(\boldsymbol{y}_i - \bar{\boldsymbol{y}})^{\textsf{T}}\) es la matriz de covarianza muestral.

     Para demostrarlo, note que la distribución previa de información unitaria satisface \(p_U(\boldsymbol{\theta}, \mathbf{\Omega}) = p_U(\boldsymbol{\theta} \mid \mathbf{\Omega}) \, p_U(\mathbf{\Omega})\), donde
     \[ \log p_U(\boldsymbol{\theta}, \mathbf{\Omega}) = \frac{1}{n} \ell(\boldsymbol{\theta}, \mathbf{\Omega} \mid \boldsymbol{y}_1,\ldots,\boldsymbol{y}_n) + \text{c}, \]
     siendo \(\ell(\boldsymbol{\theta}, \mathbf{\Omega} \mid \boldsymbol{y}_1,\ldots,\boldsymbol{y}_n)\) la log-verosimilitud reparametrizada en términos de \(\mathbf{\Omega}\), y \(\text{c}\) una constante. Además, recuerde que \((\boldsymbol{y}_i - \boldsymbol{\theta}) = (\boldsymbol{y}_i - \bar{\boldsymbol{y}}) + (\bar{\boldsymbol{y}} - \boldsymbol{\theta})\), y que una expresión de la forma \(\sum_{i=1}^n \boldsymbol{a}_i^{\textsf{T}} \mathbf{B} \boldsymbol{a}_i\) se puede escribir como \(\textsf{tr}(\mathbf{B} \mathbf{A})\), donde \(\mathbf{A} = \sum_{i=1}^n \boldsymbol{a}_i \boldsymbol{a}_i^{\textsf{T}}\).

  2. Sea \(p_U(\mathbf{\Sigma})\) la densidad inversa-Wishart inducida por \(p_U(\mathbf{\Omega})\). Obtenga una densidad proporcional a
     \[ p_U(\boldsymbol{\theta}, \mathbf{\Sigma} \mid \boldsymbol{y}_1, \dots, \boldsymbol{y}_n) \propto p(\boldsymbol{y}_1, \dots, \boldsymbol{y}_n \mid \boldsymbol{\theta}, \mathbf{\Sigma}) \, p_U(\boldsymbol{\theta} \mid \mathbf{\Sigma}) \, p_U(\mathbf{\Sigma}) . \]

• Los archivos bluecrab.dat y orangecrab.dat contienen mediciones en milímetros de la profundidad del cuerpo (\(y_1\)) y el ancho posterior (\(y_2\)) de 50 cangrejos machos de cada una de las dos especies: azul y naranja. Se modelarán estos datos utilizando una distribución Normal bivariada.

  1. Para cada una de las dos especies, obtenga las distribuciones posteriores del vector de medias poblacional \(\boldsymbol{\theta}\) y la matriz de covarianza \(\mathbf{\Sigma}\) utilizando distribuciones previas semiconjugadas. Fije \(\boldsymbol{\mu}_0\) como la media muestral de los datos, \(\mathbf{\Lambda}_0\) y \(\mathbf{S}_0\) como la matriz de covarianza muestral, y establezca \(\nu_0 = 4\). Genere 10,000 muestras posteriores de \(\boldsymbol{\theta}\) y \(\mathbf{\Sigma}\). Realice un análisis exhaustivo de convergencia.

     Cabe notar que esta elección de distribución previa centra los parámetros alrededor de estimaciones empíricas obtenidas a partir de los datos observados, similar a la distribución previa de información unitaria. Sin embargo, dado que se deriva de los datos observados, no se puede considerar una distribución previa genuina, aunque sí como una distribución previa que refleja una información inicial débil pero no sesgada.

  2. Represente gráficamente la posterior de \(\boldsymbol{\theta} = (\theta_{\text{azul}}, \theta_{\text{naranja}})\) para cada grupo y compárelos. Analice y describa cualquier diferencia en el tamaño entre las dos especies.

  3. A partir de cada matriz de covarianza obtenida mediante el muestreador de Gibbs, calcule el coeficiente de correlación correspondiente. Con estos valores, represente las densidades posteriores de las correlaciones \(\rho_{\text{azul}}\) y \(\rho_{\text{naranja}}\) para ambos grupos. Evalúe las diferencias entre las dos especies comparando estas distribuciones posteriores. En particular, estime la probabilidad \(\textsf{Pr}(\rho_{\text{azul}} < \rho_{\text{naranja}} \mid \boldsymbol{y}_{\text{azul}}, \boldsymbol{y}_{\text{naranja}})\). Interprete los resultados y analice qué sugieren sobre las diferencias entre ambas poblaciones.

• El archivo agehw.dat contiene información sobre las edades de 100 parejas casadas seleccionadas de una muestra representativa de la población de EE.UU.

  1. Antes de examinar los datos, utilice su estado de información para formular una distribución previa semiconjugada para \(\boldsymbol{\theta} = (\theta_H, \theta_W)^{\textsf{T}}\) y \(\mathbf{\Sigma}\), donde \(\theta_H\) y \(\theta_W\) representan las edades promedio de los esposos y esposas, respectivamente, y \(\mathbf{\Sigma}\) es la matriz de covarianza.

  2. Genere un conjunto de datos simulados de tamaño \(n = 100\) a partir de la distribución predictiva previa, muestreando primero \((\boldsymbol{\theta}, \mathbf{\Sigma})\) de la distribución previa y luego simulando \(\boldsymbol{y}_1, \dots, \boldsymbol{y}_n \sim \textsf{N}(\boldsymbol{\theta}, \mathbf{\Sigma})\) de manera independiente. Repita este procedimiento varias veces, genere diagramas de dispersión bivariados para cada conjunto de datos simulado y verifique si representan razonablemente sus creencias previas sobre cómo debería lucir un conjunto de datos real. Si los datos simulados no reflejan adecuadamente sus creencias, regrese al apartado anterior y reformule una nueva distribución previa. Finalmente, reporte la distribución previa elegida y presente diagramas de dispersión de al menos tres conjuntos de datos generados a partir de la distribución predictiva previa.

  3. Utilizando la distribución previa y los 100 valores del conjunto de datos, obtenga una aproximación MCMC a la distribución posterior \(p(\boldsymbol{\theta}, \mathbf{\Sigma} \mid \boldsymbol{y}_1, \dots, \boldsymbol{y}_{100})\). Represente gráficamente la distribución conjunta posterior de \(\theta_H\) y \(\theta_W\), así como la densidad marginal posterior del coeficiente de correlación entre \(y_H\) y \(y_W\), que representan las edades de los esposos y esposas, respectivamente. Además, calcule los intervalos de credibilidad posteriores del 95% para \(\theta_H\), \(\theta_W\) y el coeficiente de correlación.

  4. Calcule los intervalos de credibilidad posteriores del 95% para \(\theta_H\), \(\theta_W\) y el coeficiente de correlación utilizando las siguientes distribuciones previas:

     • La previa de Jeffreys.
     • La previa de información unitaria.
     • Una previa difusa semiconjugada, con \(\boldsymbol{\mu}_0 = \mathbf{0}\), \(\mathbf{\Lambda}_0 = 10^5 \times \mathbf{I}\), \(\mathbf{S}_0 = 1000 \times \mathbf{I}\) y \(\nu_0 = 3\).

  5. Compare los intervalos de credibilidad obtenidos en el inciso d. con aquellos calculados en el inciso c. Analice si la información previa utilizada resulta útil para la estimación de \(\boldsymbol{\theta}\) y \(\mathbf{\Sigma}\), o si alguna de las alternativas en d. es más adecuada. Además, considere el impacto de un tamaño muestral mucho menor, por ejemplo, \(n = 25\), y discuta cómo esto afectaría la influencia de la distribución previa en la estimación.

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.