Resumen Técnico: Modelos Lineales Bayesianos (BLM)

1. Planteamiento del Problema

El artículo aborda la necesidad de una presentación general, compacta y accesible de los resultados algebraicos fundamentales para el análisis bayesiano de modelos lineales. Aunque las derivaciones de la distribución posterior conjugada (Normal-Gamma) existen en numerosos libros de texto, a menudo están dispersas o presentadas de manera no unificada.

El problema central es proporcionar un marco matemático unificado que cubra:

• Modelos lineales con observaciones correlacionadas (donde la matriz de covarianza $\Sigma$ no es la identidad).
• Predicción y selección de modelos.
• Casos particulares como regresión estándar, Procesos Gaussianos (GP) y Modelos Lineales Dinámicos Bayesianos (DLM).
• La implementación eficiente de estos cálculos numéricamente, evitando la inversión directa de matrices grandes cuando es posible.

2. Metodología

El autor desarrolla una derivación algebraica rigurosa paso a paso para obtener la distribución posterior conjunta de los parámetros del modelo ($\theta$) y la precisión ($\lambda$).

2.1. El Modelo y los Prioris

Se define el modelo lineal general:
$$Y = X\theta + \epsilon; \quad \epsilon \sim N_n(0, \lambda^{-1}\Sigma)$$
Donde:

• $Y$: Vector de respuestas ($n \times 1$).
• $X$: Matriz de diseño ($n \times p$).
• $\theta$: Vector de parámetros desconocidos ($p \times 1$).
• $\Sigma$: Matriz de covarianza conocida (o su inversa $A = \Sigma^{-1}$).
• $\lambda$: Parámetro de precisión desconocido.

Se asume un prior conjugado Normal-Gamma:

• $\theta | \lambda \sim N_p(\theta_0, \lambda^{-1}A_0^{-1})$
• $\lambda \sim Ga(\alpha_0, \beta_0)$

2.2. Derivación de la Posterior

El núcleo del trabajo es la expansión de la forma cuadrática en el exponente de la función de verosimilitud combinada con el prior. Utilizando identidades matriciales clave (Lema de Sherman-Woodbury-Morrison y Lema del Determinante de Matriz), el autor demuestra que la posterior conjunta $\pi(\theta, \lambda | y)$ también sigue una distribución Normal-Gamma con parámetros actualizados:

• Precisión posterior de $\theta$: $A_n = A_0 + X'AX$
• Media posterior de $\theta$: $\theta_n = A_n^{-1}(A_0\theta_0 + X'Ay)$
• Parámetros de $\lambda$: $\alpha_n = \alpha_0 + n/2$ y una nueva tasa $\beta_n$.

El autor presenta tres expresiones equivalentes para $\beta_n$, destacando una que separa la suma de cuadrados de los residuos ponderados y la penalización por la discrepancia entre la media a priori y la posterior:
$$\beta_n = \beta_0 + \frac{1}{2} \left[ (y - X\theta_n)'A(y - X\theta_n) + (\theta_0 - \theta_n)'A_0(\theta_0 - \theta_n) \right]$$

2.3. Distribuciones Predictivas y Marginales

• Marginal de $\lambda$: Sigue una distribución Gamma $\text{Ga}(\alpha_n, \beta_n)$.
• Marginal de $\theta$: Sigue una distribución $t$ multivariada.
• Distribución Predictiva: Para nuevas observaciones $Z$ con covariables $X_m$, se deriva una distribución predictiva $t$ multivariada, incorporando la incertidumbre tanto en $\theta$ como en $\lambda$.

2.4. Eficiencia Computacional

Se proponen estrategias para optimizar el cálculo cuando $n \gg p$:

• Uso de la descomposición de Cholesky de $\Sigma = LL'$.
• Transformación de la matriz de diseño $X$ para evitar calcular $X'AX$ directamente en cada paso.
• Aplicación del Lema del Determinante de Matriz para calcular la evidencia del modelo (constante de normalización) sin invertir matrices de tamaño $n \times n$, reduciendo la complejidad computacional.

3. Contribuciones Clave

1. Unificación General: Proporciona un conjunto de fórmulas cerradas que son válidas tanto para datos independientes ($\Sigma = I$) como correlacionados ($\Sigma \neq I$), abarcando series temporales y estadística espacial.
2. Derivación Explícita de $\beta_n$: Ofrece una forma de calcular el parámetro de tasa posterior que es numéricamente estable y conceptualmente clara (suma de residuos + penalización del prior).
3. Marco para Selección de Modelos: Define la "evidencia del modelo" (marginal likelihood) necesaria para comparar modelos con diferentes matrices de diseño $X_i$ o priores, utilizando la constante de normalización derivada.
4. Implementación de Software: Presenta una clase en Python (BLM) disponible en GitHub que implementa todas las fórmulas derivadas, facilitando la aplicación práctica.
5. Cobertura de Casos Especiales: Demuestra cómo la regresión lineal estándar, los Procesos Gaussianos (GP) y los Modelos Lineales Dinámicos (DLM) son casos particulares de este marco general.

4. Resultados y Ejemplos

El autor ilustra la metodología con dos ejemplos principales:

• Regresión Polinómica: Se simulan datos de una función sinusoidal con ruido. Se ajustan modelos de regresión polinómica de grados $p=1$ a $6$. Se calculan las distribuciones marginales posteriores de los parámetros (distribuciones $t$) y las probabilidades posteriores de cada modelo, mostrando cómo el método selecciona el grado adecuado.
• Procesos Gaussianos (GP): Se aplica el marco a la interpolación de un proceso estocástico. Se demuestra cómo la predicción bayesiana en un GP es un caso particular de la distribución predictiva derivada, donde la matriz de covarianza $\Sigma$ se construye mediante una función de kernel. Se visualiza cómo el predictor actúa como interpolador (varianza cero en los puntos observados) o suavizador (si se incluye un "ruido de nugget").

5. Significado e Impacto

Este trabajo es significativo por varias razones:

• Accesibilidad: Llena un vacío al ofrecer una referencia "compacta" que evita la dispersión de resultados en múltiples textos especializados.
• Versatilidad: Al incluir observaciones correlacionadas, el marco es directamente aplicable a problemas complejos en estadística espacial y series temporales, áreas donde la independencia de las observaciones es una suposición rara.
• Eficiencia: Las notas sobre implementación numérica (uso de Cholesky y lemas matriciales) son cruciales para la viabilidad de estos modelos en conjuntos de datos grandes, donde la inversión directa de matrices sería prohibitiva.
• Reproducibilidad: La provisión de código en Python permite a investigadores y practicantes aplicar inmediatamente estas derivaciones teóricas a sus propios datos, facilitando la adopción de métodos bayesianos rigurosos.

En conclusión, el artículo establece un estándar algebraico y computacional para el análisis bayesiano de modelos lineales, demostrando que una formulación general puede simplificar y unificar el tratamiento de una amplia gama de problemas estadísticos modernos.