Tarea 6 - Regresión Lineal Multivariada Bayesiana

Alumno: Salvador Malfavón Sánchez | Profesor: Harim García Lamont | Materia: Estadística Bayesiana | Semestre: 2026-2

Consideramos el modelo de regresión lineal multivariada:

\[\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}, \qquad \boldsymbol{\varepsilon} \sim \mathcal{N}(\mathbf{0},\, \sigma^2 \mathbf{I}_n)\]

donde:

Símbolo	Dimensión	Descripción
\(\mathbf{y}\)	\(n \times 1\)	Vector de respuestas observadas
\(\mathbf{X}\)	\(n \times p\)	Matriz de diseño (incluye columna de unos para el intercepto)
\(\boldsymbol{\beta}\)	\(p \times 1\)	Vector de coeficientes desconocidos
\(\sigma^2\)	escalar	Varianza del error (desconocida)

El objetivo es obtener la distribución a posteriori de \((\boldsymbol{\beta}, \sigma^2)\) dado \(\mathbf{y}\), usando la g-prior de Zellner como distribución a priori.

1. Verosimilitud

Dado que las observaciones son independientes con \(y_i \mid \boldsymbol{\beta}, \sigma^2 \sim \mathcal{N}(\mathbf{x}_i^\top\boldsymbol{\beta},\, \sigma^2)\), la función de verosimilitud es:

\[p(\mathbf{y} \mid \boldsymbol{\beta}, \sigma^2) = \prod_{i=1}^{n} \frac{1}{\sqrt{2\pi\sigma^2}} \exp\!\left(-\frac{(y_i - \mathbf{x}_i^\top\boldsymbol{\beta})^2}{2\sigma^2}\right)\]

Agrupando en forma matricial:

\[\boxed{p(\mathbf{y} \mid \boldsymbol{\beta}, \sigma^2) = (2\pi\sigma^2)^{-n/2} \exp\!\left(-\frac{1}{2\sigma^2}\,\|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2\right)}\]

donde \(\|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2 = (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^\top(\mathbf{y} - \mathbf{X}\boldsymbol{\beta})\) es la suma de residuos al cuadrado.

Interpretación: La verosimilitud concentra su masa alrededor de los valores \(\boldsymbol{\beta}\) que minimizan la distancia cuadrática entre \(\mathbf{y}\) y \(\mathbf{X}\boldsymbol{\beta}\), es decir, alrededor del estimador de mínimos cuadrados ordinarios (MCO):

\[\hat{\boldsymbol{\beta}}_{\text{MCO}} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}\]

2. Posterior Conjunta

2.1 Distribución a priori: G-Prior de Zellner

La g-prior es una distribución normal para \(\boldsymbol{\beta}\) condicionada a \(\sigma^2\):

\[\boldsymbol{\beta} \mid \sigma^2,\, g \;\sim\; \mathcal{N}\!\left(\mathbf{0},\; g\,\sigma^2\,(\mathbf{X}^\top\mathbf{X})^{-1}\right)\]

\[p(\boldsymbol{\beta} \mid \sigma^2) = (2\pi)^{-p/2}\,\left|g\sigma^2(\mathbf{X}^\top\mathbf{X})^{-1}\right|^{-1/2} \exp\!\left(-\frac{1}{2g\sigma^2}\,\boldsymbol{\beta}^\top(\mathbf{X}^\top\mathbf{X})\boldsymbol{\beta}\right)\]

El hiperparámetro \(g > 0\) controla la difusidad del prior: - \(g\) grande → prior difusa (más peso a los datos, resultado cercano a MCO). - \(g\) pequeño → prior concentrada en \(\boldsymbol{\beta} = \mathbf{0}\) (contracción fuerte).

Propiedad clave: La covarianza a priori es proporcional a la covarianza del estimador MCO, \(\text{Var}(\hat{\boldsymbol{\beta}}_{\text{MCO}}) = \sigma^2(\mathbf{X}^\top\mathbf{X})^{-1}\). Esto hace que la g-prior sea invariante a la escala de las columnas de \(\mathbf{X}\) y que las inferencias sean coherentes.

Para \(\sigma^2\) usamos la prior no informativa de Jeffreys:

\[p(\sigma^2) \propto \frac{1}{\sigma^2}, \qquad \sigma^2 > 0\]

La distribución conjunta a priori es entonces:

\[p(\boldsymbol{\beta}, \sigma^2) \propto (\sigma^2)^{-p/2 - 1} \exp\!\left(-\frac{1}{2g\sigma^2}\,\boldsymbol{\beta}^\top(\mathbf{X}^\top\mathbf{X})\boldsymbol{\beta}\right)\]

2.2 Derivación de la posterior conjunta

Por el teorema de Bayes:

\[p(\boldsymbol{\beta}, \sigma^2 \mid \mathbf{y}) \;\propto\; p(\mathbf{y} \mid \boldsymbol{\beta}, \sigma^2)\cdot p(\boldsymbol{\beta} \mid \sigma^2)\cdot p(\sigma^2)\]

Combinando los términos:

\[p(\boldsymbol{\beta}, \sigma^2 \mid \mathbf{y}) \;\propto\; (\sigma^2)^{-(n+p)/2\,-\,1} \exp\!\left(-\frac{1}{2\sigma^2}\,Q(\boldsymbol{\beta})\right)\]

donde definimos la forma cuadrática:

\[Q(\boldsymbol{\beta}) \;=\; \|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2 + \frac{1}{g}\,\boldsymbol{\beta}^\top(\mathbf{X}^\top\mathbf{X})\boldsymbol{\beta}\]

Completación del cuadrado en \(\boldsymbol{\beta}\): Expandiendo \(Q(\boldsymbol{\beta})\):

\[Q(\boldsymbol{\beta}) = \mathbf{y}^\top\mathbf{y} - 2\boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{y} + \boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{X}\boldsymbol{\beta} + \frac{1}{g}\boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{X}\boldsymbol{\beta}\]

\[= \mathbf{y}^\top\mathbf{y} - 2\boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{y} + \frac{g+1}{g}\,\boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{X}\boldsymbol{\beta}\]

Definimos la matriz de precisión posterior (dividida por \(\sigma^2\)) y completamos el cuadrado:

\[Q(\boldsymbol{\beta}) = (\boldsymbol{\beta} - \boldsymbol{\beta}_g)^\top \mathbf{V}_g^{-1}(\boldsymbol{\beta} - \boldsymbol{\beta}_g) + S_g\]

donde los parámetros posteriores son:

\[\boldsymbol{\beta}_g = \frac{g}{g+1}\,\hat{\boldsymbol{\beta}}_{\text{MCO}}\]

\[\mathbf{V}_g = \frac{g}{g+1}\,(\mathbf{X}^\top\mathbf{X})^{-1}\]

\[S_g = \mathbf{y}^\top\mathbf{y} - \frac{g}{g+1}\,\hat{\boldsymbol{\beta}}_{\text{MCO}}^\top\mathbf{X}^\top\mathbf{y}\]

Verificación de la completación: La inversa de \(\mathbf{V}_g\) es \(\mathbf{V}_g^{-1} = \frac{g+1}{g}(\mathbf{X}^\top\mathbf{X})\). Entonces \(\mathbf{V}_g^{-1}\boldsymbol{\beta}_g = \frac{g+1}{g}(\mathbf{X}^\top\mathbf{X})\cdot\frac{g}{g+1}\hat{\boldsymbol{\beta}}_{\text{MCO}} = (\mathbf{X}^\top\mathbf{X})\hat{\boldsymbol{\beta}}_{\text{MCO}} = \mathbf{X}^\top\mathbf{y}\). ✓

2.3 Distribuciones a posteriori marginales

Sustituyendo la completación del cuadrado en la posterior conjunta:

\[p(\boldsymbol{\beta}, \sigma^2 \mid \mathbf{y}) \;\propto\; (\sigma^2)^{-(n+p)/2\,-\,1} \exp\!\left(-\frac{(\boldsymbol{\beta} - \boldsymbol{\beta}_g)^\top \mathbf{V}_g^{-1}(\boldsymbol{\beta} - \boldsymbol{\beta}_g) + S_g}{2\sigma^2}\right)\]

Posterior condicional de \(\boldsymbol{\beta}\) dado \(\sigma^2\): Identificando el kernel de una Normal multivariada en \(\boldsymbol{\beta}\):

\[\boxed{\boldsymbol{\beta} \mid \sigma^2, \mathbf{y} \;\sim\; \mathcal{N}\!\left(\boldsymbol{\beta}_g,\; \sigma^2\,\mathbf{V}_g\right)}\]

Posterior marginal de \(\sigma^2\): Integrando \(\boldsymbol{\beta}\) (integral gaussiana multivariada, que aporta un factor \((\sigma^2)^{p/2}\)):

\[p(\sigma^2 \mid \mathbf{y}) \;\propto\; (\sigma^2)^{-n/2\,-\,1} \exp\!\left(-\frac{S_g}{2\sigma^2}\right)\]

\[\boxed{\sigma^2 \mid \mathbf{y} \;\sim\; \text{Inv-Gamma}\!\left(\frac{n}{2},\; \frac{S_g}{2}\right)}\]

Posterior marginal de \(\boldsymbol{\beta}\): Integrando \(\sigma^2\) de la posterior conjunta, se obtiene una distribución \(t\) de Student multivariada:

\[\boxed{\boldsymbol{\beta} \mid \mathbf{y} \;\sim\; t_n\!\left(\boldsymbol{\beta}_g,\; \frac{S_g}{n}\,\mathbf{V}_g\right)}\]

con \(n\) grados de libertad, vector de localización \(\boldsymbol{\beta}_g\) y matriz de escala \(\frac{S_g}{n}\mathbf{V}_g\).

La posterior conjunta \((\boldsymbol{\beta}, \sigma^2) \mid \mathbf{y}\) pertenece a la familia Normal-Gamma Inversa (NIG), que es la familia conjugada para este modelo.

3. Predictiva

La distribución predictiva a posteriori describe la distribución de una nueva observación \(\tilde{y}^*\) en un nuevo punto \(\mathbf{x}^*\), incorporando toda la incertidumbre sobre \((\boldsymbol{\beta}, \sigma^2)\):

\[p(\tilde{y}^* \mid \mathbf{y}) = \iint p(\tilde{y}^* \mid \boldsymbol{\beta}, \sigma^2)\; p(\boldsymbol{\beta}, \sigma^2 \mid \mathbf{y})\; d\boldsymbol{\beta}\, d\sigma^2\]

Paso 1 — Integrar sobre \(\boldsymbol{\beta}\): Como \(\tilde{y}^* \mid \boldsymbol{\beta}, \sigma^2 \sim \mathcal{N}({\mathbf{x}^*}^\top\boldsymbol{\beta}, \sigma^2)\) y \(\boldsymbol{\beta} \mid \sigma^2, \mathbf{y} \sim \mathcal{N}(\boldsymbol{\beta}_g, \sigma^2\mathbf{V}_g)\), la marginalización produce:

\[\tilde{y}^* \mid \sigma^2, \mathbf{y} \;\sim\; \mathcal{N}\!\left({\mathbf{x}^*}^\top\boldsymbol{\beta}_g,\; \sigma^2\underbrace{\left(1 + {\mathbf{x}^*}^\top\mathbf{V}_g\mathbf{x}^*\right)}_{c^*}\right)\]

El término \(1\) corresponde al ruido de predicción \(\varepsilon^*\), y \({\mathbf{x}^*}^\top\mathbf{V}_g\mathbf{x}^*\) a la incertidumbre sobre \(\boldsymbol{\beta}\).

Paso 2 — Integrar sobre \(\sigma^2\): Como \(\sigma^2 \mid \mathbf{y} \sim \text{Inv-Gamma}(n/2, S_g/2)\), la mezcla de normales con varianza Inv-Gamma produce una distribución \(t\):

\[\boxed{\tilde{y}^* \mid \mathbf{y} \;\sim\; t_n\!\left({\mathbf{x}^*}^\top\boldsymbol{\beta}_g,\; \frac{S_g}{n}\left(1 + {\mathbf{x}^*}^\top\mathbf{V}_g\mathbf{x}^*\right)\right)}\]

con \(n\) grados de libertad.

Interpretación de los componentes:

\({\mathbf{x}^*}^\top\boldsymbol{\beta}_g\): predicción puntual (media a posteriori contraída hacia cero).
\(\frac{S_g}{n}\): estimación posterior de \(\sigma^2\) (varianza del ruido).
\({\mathbf{x}^*}^\top\mathbf{V}_g\mathbf{x}^*\): incertidumbre en la estimación de \(\boldsymbol{\beta}\) (análogo al error estándar de predicción en MCO).
Grados de libertad \(n\): colas más pesadas que la normal, reflejando incertidumbre paramétrica.

La varianza de la predictiva es:

\[\text{Var}(\tilde{y}^* \mid \mathbf{y}) = \frac{n}{n-2}\cdot\frac{S_g}{n}\left(1 + {\mathbf{x}^*}^\top\mathbf{V}_g\mathbf{x}^*\right) = \frac{S_g}{n-2}\left(1 + {\mathbf{x}^*}^\top\mathbf{V}_g\mathbf{x}^*\right), \quad n > 2\]

4. Estimadores Bayesianos

Un estimador Bayesiano de un parámetro \(\theta\) es el resumen puntual de su distribución a posteriori que minimiza una función de pérdida \(L(\hat{\theta}, \theta)\).

4.1 Estimadores para \(\boldsymbol{\beta}\)

La posterior marginal de \(\boldsymbol{\beta}\) es una \(t\) multivariada simétrica y unimodal centrada en \(\boldsymbol{\beta}_g\). Por lo tanto, media = mediana = moda:

\[\boxed{\hat{\boldsymbol{\beta}}_{\text{Bayes}} = \mathbb{E}[\boldsymbol{\beta} \mid \mathbf{y}] = \boldsymbol{\beta}_g = \frac{g}{g+1}\,\hat{\boldsymbol{\beta}}_{\text{MCO}}}\]

Este estimador tiene tres propiedades notables:

Contracción (shrinkage): El factor \(\frac{g}{g+1} < 1\) encoge las estimaciones hacia el vector de medias a priori \(\mathbf{0}\). Esto reduce varianza al costo de introducir un pequeño sesgo.
Casos extremos:
- \(g \to \infty\): \(\hat{\boldsymbol{\beta}}_{\text{Bayes}} \to \hat{\boldsymbol{\beta}}_{\text{MCO}}\) (prior no informativa, domina la verosimilitud).
- \(g \to 0\): \(\hat{\boldsymbol{\beta}}_{\text{Bayes}} \to \mathbf{0}\) (prior domina, colapsa al valor a priori).
Conexión con regularización: Para \(g = (n-1)(p+1)/\lambda\) se recupera la regresión ridge con parámetro \(\lambda\), lo que muestra que la g-prior tiene una interpretación como regularizador \(L_2\).

Funciones de pérdida asociadas:

Función de pérdida	Estimador óptimo
Cuadrática: \(L = \\|\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}\\|^2\)	Media posterior: \(\boldsymbol{\beta}_g\)
Valor absoluto: \(L = \\|\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}\\|_1\)	Mediana posterior: \(\boldsymbol{\beta}_g\) (coincide)
0-1 (error de decisión): \(L = \mathbf{1}[\hat{\boldsymbol{\beta}} \neq \boldsymbol{\beta}]\)	Moda posterior: \(\boldsymbol{\beta}_g\) (coincide)

4.2 Estimadores para \(\sigma^2\)

La posterior es \(\sigma^2 \mid \mathbf{y} \sim \text{Inv-Gamma}(n/2,\, S_g/2)\). Sus resúmenes son:

\[\mathbb{E}[\sigma^2 \mid \mathbf{y}] = \frac{S_g/2}{n/2 - 1} = \frac{S_g}{n - 2}, \quad n > 2 \qquad \text{(media posterior)}\]

\[\text{Mo}[\sigma^2 \mid \mathbf{y}] = \frac{S_g/2}{n/2 + 1} = \frac{S_g}{n + 2} \qquad \text{(moda posterior = MAP)}\]

4.3 Elección del hiperparámetro \(g\)

La elección de \(g\) es crucial. Algunas opciones comunes son:

Nombre	Valor	Descripción
Unit Information Prior	\(g = n\)	El prior contiene la misma información que una sola observación. Opción estándar y bien calibrada.
RIC (Kass-Wasserman)	\(g = p^2\)	Criterio de inflación de riesgo; favorece modelos parsimoniosos.
BRIC	\(g = \max(n, p^2)\)	Combina ambos criterios.
Empírico Bayes	\(\hat{g} = \arg\max_g p(\mathbf{y} \mid g)\)	Maximiza la verosimilitud marginal; adapta \(g\) a los datos.

5. Implementación del modelo en Stan con datos

5.1 Código Stan

El siguiente modelo Stan implementa la regresión lineal con g-prior:

data {
  int<lower=1> n;                   // número de observaciones
  int<lower=1> p;                   // número de columnas de X (con intercepto)
  vector[n] y;                      // variable respuesta
  matrix[n, p] X;                   // matriz de diseño
  real<lower=0> g;                  // hiperparámetro de la g-prior
}

transformed data {
  matrix[p, p] XtX     = X' * X;           // X'X
  matrix[p, p] XtX_inv = inverse_spd(XtX); // (X'X)^{-1}
  vector[p]    mu0     = rep_vector(0.0, p); // media de la g-prior
}

parameters {
  vector[p]        beta;   // coeficientes de regresión
  real<lower=0> sigma2;    // varianza del error
}

model {
  // G-prior de Zellner: beta | sigma2 ~ N(0, g * sigma2 * (X'X)^{-1})
  beta ~ multi_normal(mu0, g * sigma2 * XtX_inv);

  // Prior de Jeffreys para sigma2: p(sigma2) ∝ 1/sigma2
  target += -log(sigma2);

  // Verosimilitud: y | beta, sigma2 ~ N(X*beta, sigma2 * I_n)
  y ~ normal(X * beta, sqrt(sigma2));
}

generated quantities {
  vector[n] y_rep;
  vector[n] mu = X * beta;

  for (i in 1:n)
    y_rep[i] = normal_rng(mu[i], sqrt(sigma2));
}

5.2 Código R: simulación y ajuste

library(rstan)
library(ggplot2)
library(bayesplot)

options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

# ── 1. Simular datos ──────────────────────────────────────────────────────────
set.seed(42)

n <- 50                           # observaciones
p <- 3                            # columnas de X (intercepto + 2 predictores)
g <- n                            # unit information prior: g = n

beta_true   <- c(1.0, 2.0, -1.5) # coeficientes verdaderos
sigma2_true <- 1.0

X <- cbind(1, matrix(rnorm(n * (p - 1)), nrow = n))  
y <- as.vector(X %*% beta_true + rnorm(n, 0, sqrt(sigma2_true)))

# ── 2. Estimadores analíticos ─────────────────────────────────────────────────
XtX   <- t(X) %*% X
beta_ols <- solve(XtX) %*% t(X) %*% y

beta_g_analitico   <- g / (g + 1) * beta_ols
S_g                <- as.numeric(t(y) %*% y - g/(g+1) * t(beta_ols) %*% XtX %*% beta_ols)
sigma2_g_analitico <- S_g / (n - 2)

cat("=== Estimadores analíticos ===\n")
cat(sprintf("β_g  (Bayes): [%.3f, %.3f, %.3f]\n",
            beta_g_analitico[1], beta_g_analitico[2], beta_g_analitico[3]))
cat(sprintf("β̂   (MCO):   [%.3f, %.3f, %.3f]\n",
            beta_ols[1], beta_ols[2], beta_ols[3]))
cat(sprintf("β   (true):  [%.3f, %.3f, %.3f]\n",
            beta_true[1], beta_true[2], beta_true[3]))
cat(sprintf("σ²_g (Bayes): %.3f  |  σ² (true): %.3f\n",
            sigma2_g_analitico, sigma2_true))

# ── 3. Ajustar el modelo Stan ─────────────────────────────────────────────────
stan_data <- list(n = n, p = p, y = y, X = X, g = g)

fit <- stan(
  file    = "modelo_g_prior.stan",   # archivo Stan
  data    = stan_data,
  chains  = 4,
  iter    = 3000,
  warmup  = 1000,
  seed    = 42
)

# ── 5. Resumen del ajuste ─────────────────────────────────────────────────────
print(fit, pars = c("beta", "sigma2"), probs = c(0.025, 0.5, 0.975))

# ── 6. Comparar medias MCMC vs. estimadores analíticos ───────────────────────
posterior_df <- as.data.frame(fit)

beta_mcmc   <- colMeans(posterior_df[, paste0("beta[", 1:p, "]")])
sigma2_mcmc <- mean(posterior_df[["sigma2"]])

cat("\n=== Comparación MCMC vs. Analítico ===\n")
cat(sprintf("β_g  (MCMC):      [%.3f, %.3f, %.3f]\n", beta_mcmc[1], beta_mcmc[2], beta_mcmc[3]))
cat(sprintf("β_g  (analítico): [%.3f, %.3f, %.3f]\n", beta_g_analitico[1], beta_g_analitico[2], beta_g_analitico[3]))
cat(sprintf("σ²   (MCMC):      %.3f\n", sigma2_mcmc))
cat(sprintf("σ²   (analítico): %.3f\n", sigma2_g_analitico))

# ── 7. Gráficas ───────────────────────────────────────────────────────────────

# 7a. Distribuciones a posteriori de cada β_j
mcmc_areas(
  posterior_df[, paste0("beta[", 1:p, "]")],
  prob = 0.95
) +
  geom_vline(xintercept = beta_true, linetype = "dashed", color = "red") +
  labs(title = "Posterior de β con g-prior (línea roja = valor verdadero)",
       x = "Valor", y = "Densidad")

# 7b. Posterior predictive check: y_rep vs. y observado
y_rep_mat <- as.matrix(fit, pars = "y_rep")
ppc_dens_overlay(y, y_rep_mat[1:200, ]) +
  labs(title = "Verificación predictiva a posteriori",
       subtitle = "y observado (negro) vs. réplicas y_rep (azul)")

# 7c. Comparación de estimadores
tabla <- data.frame(
  Parámetro  = c("β₀", "β₁", "β₂"),
  Verdadero  = beta_true,
  MCO        = round(beta_ols, 3),
  Bayesiano  = round(beta_g_analitico, 3),
  MCMC_media = round(beta_mcmc, 3)
)
print(tabla, row.names = FALSE)

5.3 Resultados

Se ejecutó el código con set.seed(42), \(n = 50\), \(p = 3\) y \(g = n = 50\). La tabla siguiente reporta los estimadores reales obtenidos:

Parámetro	Verdadero	MCO \(\hat{\boldsymbol{\beta}}\)	Bayesiano \(\boldsymbol{\beta}_g\)	Contracción
\(\beta_0\)	1.0000	0.8497	0.8331	\(\times\,0.9804\)
\(\beta_1\)	2.0000	2.0882	2.0472	\(\times\,0.9804\)
\(\beta_2\)	−1.5000	−1.4785	−1.4495	\(\times\,0.9804\)

Parámetro	Verdadero	\(S^2\) MCO	Media posterior	Moda posterior
\(\sigma^2\)	1.0000	0.8869	1.0610	0.9794

El factor de contracción \(g/(g+1) = 50/51 \approx 0.9804\) es pequeño porque \(g = n\) es grande, de modo que el estimador Bayesiano apenas se aleja del MCO. La media posterior de \(\sigma^2 = S_g/(n-2) = 1.061\) cubre el valor verdadero \(\sigma^2 = 1\), y la moda \(S_g/(n+2) = 0.979\) es más conservadora.

Figura 1 — Distribución a posteriori marginal de cada \(\beta_j\)

La región sombreada cubre el 95% de probabilidad posterior. La línea roja punteada marca el valor verdadero y la línea azul el estimador MCO, ambos caen dentro del IC 95% en los tres coeficientes.

Figura 2 — Verificación predictiva a posteriori (PPC)

La distribución de \(\mathbf{y}\) observado (negro) sigue la forma de las 150 réplicas \(\mathbf{y}^{\text{rep}}\) simuladas (azul), lo que indica que el modelo captura bien la estructura de los datos.

Figura 3 — Contracción bayesiana

El gráfico ilustra el efecto de contracción: los estimadores bayesianos (triángulos rojos) se ubican sistemáticamente entre el MCO (círculos azules) y el origen \(\mathbf{0}\), en la proporción exacta \(g/(g+1)\).

Referencias

Zellner, A. (1986). On assessing prior distributions and Bayesian regression analysis with g-prior distributions. En P. K. Goel y A. Zellner (Eds.), Bayesian Inference and Decision Techniques (pp. 233–243). North-Holland.
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., y Rubin, D. B. (2013). Bayesian Data Analysis (3.ª ed.). CRC Press. (Cap. 14–15)
Liang, F., Paulo, R., Molina, G., Clyde, M. A., y Berger, J. O. (2008). Mixtures of g priors for Bayesian variable selection. Journal of the American Statistical Association, 103(481), 410–423.
Stan Development Team. (2023). Stan Modeling Language Users Guide and Reference Manual. Versión 2.32. https://mc-stan.org