Tabla de contenido

Librarías

library(cmdstanr) #  cmdstan_model

library(posterior) # as_draws_df, mcmc_hist

library(ggpmisc)  # stat_poly_line
library(ggplot2)
library(patchwork)  
library(tinytex)

library(nortest) # muestra las funciones del paquete
library(coda)


library(flextable) # flextable
library(tidyverse) # %

Modelo frecuentista

El modelo de regresión lineal simple describe la relación entre una variable dependiente continua \(y\) y una única variable independiente o predictora \(x\).

La relación poblacional se expresa mediante la ecuación:

\[y_i = \beta_0 + \beta_1 x_i + \varepsilon_i, \quad i = 1, 2, \dots, n\]

Donde: * \(\beta_0\) es el intercepto de la línea de regresión. * \(\beta_1\) es la pendiente de la línea de regresión. * \(\varepsilon_i\) es el término de error aleatorio o perturbación para la observación \(i\).

Supuestos del modelo

  1. Linealidad: La relación entre las variables es lineal en los parámetros.

  2. Exogeneidad: El valor esperado de los errores es cero, \(E[\varepsilon_i] = 0\).

  3. Homocedasticidad: La varianza de los errores es constante, \(\text{Var}(\varepsilon_i) = \sigma^2\).

  4. No autocorrelación: Los errores son independientes entre sí, \(\text{Cov}(\varepsilon_i, \varepsilon_j) = 0\) para todo \(i \neq j\).

  5. Normalidad: Los errores siguen una distribución normal, \(\varepsilon_i \sim \mathcal{N}(0, \sigma^2)\).

Utilizando el método de Mínimos Cuadrados Ordinarios, los estimadores muestrales \(\hat{\beta}_0\) y \(\hat{\beta}_1\), están dados por

\[\hat{\beta}_1 = \frac{\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})}{\sum_{i=1}^n (x_i - \bar{x})^2} = \frac{S_{xy}}{S_{xx}}\]

\[\hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}\]

Dado que las respuestas \(y_i\) son combinaciones lineales de los errores normales \(\varepsilon_i\), los estimadores de los parámetros también se distribuyen normalmente:

Distribución de la pendiente: \[\hat{\beta}_1 \sim \mathcal{N}\left(\beta_1, \sigma_{\hat{\beta}_1}^2\right)\] donde \[\sigma_{\hat{\beta}_1}^2 = \frac{\sigma^2}{\sum_{i=1}^n (x_i - \bar{x})^2}\]

Distribución del intercepto: \[\hat{\beta}_0 \sim \mathcal{N}\left(\beta_0, \sigma_{\hat{\beta}_0}^2\right)\] donde \[\sigma_{\hat{\beta}_0}^2 = \sigma^2 \left[ \frac{1}{n} + \frac{\bar{x}^2}{\sum_{i=1}^n (x_i - \bar{x})^2} \right]\] Dado que la varianza poblacional \(\sigma^2\) es usualmente desconocida, se estima de forma insesgada a través de la varianza residual muestral (\(s^2\)):

\[s^2 = \frac{\sum_{i=1}^n (y_i - \hat{y}_i)^2}{n - 2} = \frac{\text{SCE}}{n - 2}\]

Donde \(\text{SCE}\) es la Suma de Cuadrados de los Errores (o Residuales), y la distribución del estimador está ligada a una distribución Chi-cuadrada:

\[\frac{(n-2)s^2}{\sigma^2} \sim \chi^2_{n-2}\]

Al sustituir la varianza poblacional desconocida \(\sigma^2\) por su estimador insesgado \(s^2\), se debe utilizar la distribución \(t\) de Student con \(n - 2\) grados de libertad en lugar de la distribución normal estándar.

Definimos los Errores Estándar estimados como:

\[\text{EE}(\hat{\beta}_1) = \frac{s}{\sqrt{\sum_{i=1}^n (x_i - \bar{x})^2}}\]

\[\text{EE}(\hat{\beta}_0) = s \sqrt{\frac{1}{n} + \frac{\bar{x}^2}{\sum_{i=1}^n (x_i - \bar{x})^2}}\]

Para un nivel de confianza del \((1 - \alpha) \times 100\%\), el valor crítico se denota como \(t_{\alpha/2, \, n-2}\). Los intervalos de confianza correspondientes son:

Intervalo de confianza para la pendiente: \[\text{IC}(1-\alpha) = \left[ \hat{\beta}_1 - t_{\alpha/2, \, n-2} \cdot \text{EE}(\hat{\beta}_1), \;\; \hat{\beta}_1 + t_{\alpha/2, \, n-2} \cdot \text{EE}(\hat{\beta}_1) \right]\]

Intervalo de confianza para el intercepto: \[\text{IC}(1-\alpha) = \left[ \hat{\beta}_0 - t_{\alpha/2, \, n-2} \cdot \text{EE}(\hat{\beta}_0), \;\; \hat{\beta}_0 + t_{\alpha/2, \, n-2} \cdot \text{EE}(\hat{\beta}_0) \right]\]

Las pruebas de hipótesis evalúan si existe una relación estadísticamente significativa entre las variables y si los coeficientes difieren de un valor teórico (generalmente cero).

Esta es la prueba más importante del modelo. Evalúa si la variable independiente \(x\) tiene un efecto lineal sobre \(y\).

Evalúa si la línea de regresión corta el eje vertical en un punto diferente de cero.

Cuando queremos usar el modelo estimado para estimar o predecir el valor de la variable dependiente ante un nuevo valor específico del predictor (\(x_0\)), debemos distinguir entre dos conceptos:

  1. Valor Esperado o Medio (\(E[y_0 \mid x_0]\)): El promedio de todos los valores de \(y\) que corresponden a \(x_0\). Se asocia a un Intervalo de Confianza para la Media.
  2. Valor Individual (\(y_0\)): El valor específico que tomará una nueva observación individual. Se asocia a un Intervalo de Predicción.

Dado que una observación individual incluye la variabilidad inherente del término de error aleatorio (\(\sigma^2\)), el intervalo de predicción siempre es más ancho que el intervalo de confianza.

El error estándar estimado para predecir una observación individual \(y_0\) en el punto \(x_0\) se define matemáticamente como:

\[\text{EE}(\text{pred}) = s \sqrt{1 + \frac{1}{n} + \frac{(x_0 - \bar{x})^2}{\sum_{i=1}^n (x_i - \bar{x})^2}}\]

Nota: el término \(+1\) dentro de la raíz, el cual representa la variabilidad no explicada por el modelo \(\sigma^2\).

Para un valor dado \(x_0\), el valor puntual predicho es \(\hat{y}_0 = \hat{\beta}_0 + \hat{\beta}_1 x_0\). Con un nivel de confianza del \((1 - \alpha) \times 100\%\), el intervalo de predicción para la nueva observación es:

\[\text{IP}(1-\alpha) = \left[\hat{y}_0 - t_{\alpha/2, \, n-2} \cdot \text{EE}(\text{pred}), \;\; \hat{y}_0 + t_{\alpha/2, \, n-2} \cdot \text{EE}(\text{pred})\right]\]

Como complemento a las pruebas \(t\), la significancia global del modelo se evalúa mediante un Análisis de Varianza (ANOVA), el cual descompone la variabilidad total de los datos (\(\text{SCT}\)) en la variabilidad explicada por la regresión (\(\text{SCR}\)) y la variabilidad de los residuos o errores (\(\text{SCE}\)):

\[\text{SCT} = \text{SCR} + \text{SCE}\]

\[\sum_{i=1}^n (y_i - \bar{y})^2 = \sum_{i=1}^n (\hat{y}_i - \bar{y})^2 + \sum_{i=1}^n (y_i - \hat{y}_i)^2\]

El estadístico de prueba global es una distribución \(F\) de Snedecor con \(1\) y \(n-2\) grados de libertad:

\[F_{\text{calc}} = \frac{\text{CMR}}{\text{CME}} = \frac{\text{SCR} / 1}{\text{SCE}/(n - 2)}\]

En una regresión lineal simple, el resultado de rechazar la hipótesis nula en la prueba \(F\) es idéntico al de la prueba \(t\) de la pendiente, cumpliéndose rigurosamente que \(F_{\text{calc}} = (t_{\text{calc}})^2\).

Determinar: ¿A cuánto ascenderá el salario anual de una persona que haya empezado a trabajar en esta empresa con 10 años de experiencia laboral?

ruta_dat<-"C:\\Documetos MGM_OS_Dell\\Documentos Maria_Dell\\Facultad_Matematicas\\1. Programa_Licenciatura en Matematicas\\Cursos de la LM\\Modelacion_bayesiana_LM\\Semestre_2_2025\\Bayesian_Statistical_Modeling_with_Stan_R_and_Python-master\\chap04\\input\\data-salary.csv" 

d <-read.csv(file=ruta_dat) # Datos

De acuerdo con los resultado, se tiene que los años de experiencia explican el ingreso anual del trabajador.

res_lm<-lm(Y~X,data=d)
summary(res_lm) 
## 
## Call:
## lm(formula = Y ~ X, data = d)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0063 -2.0278  0.2412  2.3603  2.9080 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 38.69687    1.33429  29.002 3.35e-13 ***
## X            0.75238    0.08261   9.107 5.26e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.614 on 13 degrees of freedom
## Multiple R-squared:  0.8645, Adjusted R-squared:  0.8541 
## F-statistic: 82.95 on 1 and 13 DF,  p-value: 5.264e-07

Note que el coeficiente de determinación tiene un valor de 0.86, lo cual indica un buen ajuste del modelo.

plot(res_lm)

hist(res_lm$residuals) # histograma de los residuales

ad.test(res_lm$residuals) # prueba de normalidad
## 
##  Anderson-Darling normality test
## 
## data:  res_lm$residuals
## A = 0.49405, p-value = 0.1824

Tarea: evaluar todos los supuestos del modelo ajustado, y generar el gráfico de los datos con el modelo ajustado

ggplot(data = d, aes(x =X, y = Y)) +
  stat_poly_line(se=FALSE)+ # Gráfica la línea
  stat_poly_eq(use_label("eq")) +
  stat_poly_eq(label.y = 0.9) +
  geom_point()+
  theme_bw()

Intervalos de confianza de los coeficientes de regresión:

confint(res_lm)
##                  2.5 %    97.5 %
## (Intercept) 35.8143149 41.579424
## X            0.5739065  0.930849

Tarea: calcular el intervalo de confianza de \(\sigma\).

Los intervalos de confianza del 95 % de los ingresos iniciales, \(a + bX_{pred}\), son calculados para los empleados con \(X_{pred}\) de los años de experiencia laboral.

X_pred<-data.frame(X=0:28)

conf_95 <-predict(res_lm, X_pred, interval="confidence", level=0.95)

data.frame(conf_95) %>% 
flextable()

fit

lwr

upr

38.69687

35.81431

41.57942

39.44925

36.71916

42.17933

40.20163

37.62067

42.78258

40.95400

38.51824

43.38977

41.70638

39.41110

44.00166

42.45876

40.29835

44.61917

43.21114

41.17886

45.24341

43.96351

42.05129

45.87574

44.71589

42.91401

46.51777

45.46827

43.76514

47.17140

46.22065

44.60255

47.83874

46.97303

45.42399

48.52206

47.72540

46.22724

49.22357

48.47778

47.01040

49.94516

49.23016

47.77222

50.68810

49.98254

48.51227

51.45280

50.73491

49.23109

52.23874

51.48729

49.93005

53.04453

52.23967

50.61111

53.86823

52.99205

51.27649

54.70760

53.74443

51.92846

55.56039

54.49680

52.56909

56.42451

55.24918

53.20026

57.29811

56.00156

53.82353

58.17959

56.75394

54.44023

59.06765

57.50631

55.05144

59.96119

58.25869

55.65807

60.85931

59.01107

56.26084

61.76130

59.76345

56.86034

62.66656

Los intervalos de confianza del 95 % de los ingresos iniciales, \(a + bX_{pred}+\epsilon\), son calculados para los empleados con \(X_{pred}\) de los años de experiencia laboral.

pred_95 <-predict(res_lm, X_pred, interval="prediction", level=0.95)

data.frame(pred_95)%>% 
flextable()

fit

lwr

upr

38.69687

32.35725

45.03649

39.44925

33.17748

45.72102

40.20163

33.99332

46.40993

40.95400

34.80464

47.10337

41.70638

35.61130

47.80146

42.45876

36.41317

48.50434

43.21114

37.21015

49.21212

43.96351

38.00211

49.92492

44.71589

38.78896

50.64282

45.46827

39.57061

51.36593

46.22065

40.34698

52.09431

46.97303

41.11801

52.82804

47.72540

41.88364

53.56717

48.47778

42.64384

54.31173

49.23016

43.39858

55.06174

49.98254

44.14786

55.81721

50.73491

44.89170

56.57813

51.48729

45.63010

57.34449

52.23967

46.36311

58.11623

52.99205

47.09079

58.89330

53.74443

47.81320

59.67565

54.49680

48.53042

60.46319

55.24918

49.24253

61.25583

56.00156

49.94965

62.05346

56.75394

50.65189

62.85598

57.50631

51.34936

63.66327

58.25869

52.04219

64.47520

59.01107

52.73050

65.29164

59.76345

53.41445

66.11244

Modelo bayesiano

Desarrollo analítico

Stan es generalmente mejor si buscas estabilidad, algoritmos MCMC ultra optimizados de fábrica y una comunidad masiva, mientras que Julia (a través de su ecosistema de modelado probabilístico, principalmente el paquete Turing.jl) es mejor si necesitas flexibilidad matemática extrema, personalización de algoritmos o integración nativa en producción sin cambiar de lenguaje.

Característica Stan (CmdStan / RStan / PyStan) Julia (Turing.jl / DynamicPPL)
Algoritmo NUTS (No-U-Turn Sampler) NUTS, HMC, MH, SMC

La distribución posterior es proporcional al producto de la función de verosimilitud más la distribución a priori:

\[ p(\theta|y) \propto p(Y|\theta)p(\theta) \] - Logaritmo de la distribución posterior (lp_): \[ \log p(Y|\theta)+\log p(\theta) \] - Intervalo de confianza bayesiano (intervalo de credibilidad): El intervalo comprendido entre \(\alpha/2\) y \((1-\alpha)/2\) de la distribución a posteriori se denomina \((1-\alpha)\%\) intervalo de confianza bayesiano.

  • Una de las características más importantes de Stan es que utiliza el No-U-Turn Sampler (NUTS) como algoritmo predeterminado para calcular la estimación. NUTS es una implementación del método de Monte Carlo hamiltoniano (HMC), que es un tipo de MCMC. El punto fuerte de NUTS es que permite realizar un muestreo eficaz incluso cuando el número de parámetros es elevado. En comparación con WinBUGS y JAGS, una iteración de muestreo con NUTS lleva más tiempo debido a la complejidad de su algoritmo, pero la autocorrelación entre iteraciones es menor. Un cálculo que con MCMC que requiere 100 000 iteraciones en WinBUGS o JAGS, en Stan solo se requieren 1 000 iteraciones.

En el contexto del enfoque bayesiano, la inferencia sobre una regresión lineal simple no arroja un único conjunto de valores óptimos, sino una distribución posterior para sus parámetros.

Si asumimos que los datos siguen el modelo clásico:

\[ y_i = \beta_0 + \beta_1 x_i + \varepsilon_i, \quad \text{donde } \varepsilon_i \sim \mathcal{N}(0, \sigma^2) \]

En notación vectorial, se tiene \[ \boldsymbol{\beta} = [\beta_0, \beta_1]^t \] En notación matricial, se tiene que \[\mathbf{y} \sim \mathcal{N}(\mathbf{X}\boldsymbol{\beta}, \sigma^2 I) \] \(\mathbf{X}\) es la matriz de diseño de tamaño \(n \times 2\) (una columna de unos para el intercepto \(\beta_0\) y una columna con la variable explicativa para la pendiente \(\beta_1\), e \(\mathbf{y}\) es el vector de respuestas de tamaño \(n \times 1\).

La función de verosimilitud está dada por
\[ p(\mathbf{y} \vert \mathbf{X}, \boldsymbol{\beta}, \sigma) \propto \exp\left( -\frac{1}{2\sigma^2} (\mathbf{y} - \mathbf{X}\beta)^t(\mathbf{y} - \mathbf{X}\beta) \right) \]

Considerando las distribuciones a priori conjugadas

\[ \boldsymbol{\beta} \sim \mathcal{N}(\boldsymbol{\mu}_0, \boldsymbol{\Sigma}_0) \]

\[ \sigma \sim \text{Inv-Gamma}(a_0, b_0) \]

Para resolver el modelo de forma conjunta cuando \(\sigma\) es desconocida, se suele estructurar de manera condicional:

\[\boldsymbol{\beta} \mid \sigma \sim \mathcal{N}(\boldsymbol{\mu}_0, \boldsymbol{\Sigma}_0)\] \[\sigma \sim \text{Inv-Gamma}(a_0, b_0)\] Por lo tanto: Distribución a priori para el vector \(\boldsymbol{\beta}\) \[ p(\boldsymbol{\beta} \mid \sigma) = (2\pi)^{-\frac{k}{2}} |\boldsymbol{\Sigma}_0|^{-\frac{1}{2}} \exp\left( -\frac{1}{2} (\boldsymbol{\beta} - \boldsymbol{\mu}_0)^T \boldsymbol{\Sigma}_0^{-1} (\boldsymbol{\beta} - \boldsymbol{\mu}_0) \right) \] Distribución a priori para \(\sigma^2\) \[ p(\sigma) = \frac{b_0^{a_0}}{\Gamma(a_0)} (\sigma^2)^{-(a_0 + 1)} \exp\left( -\frac{b_0}{\sigma^2} \right) \]

El producto de las función de verosimilitud y las distribuciones a priori, se tiene la distribución posterior conjunta (Normal-Inversa-Gamma) de los coeficientes \(\boldsymbol{\beta}\) y la desviación estándar \(\sigma\), dados los datos \((\mathbf{y}, \mathbf{X})\):

\[p(\boldsymbol{\beta}, \sigma \mid \mathbf{y}, \mathbf{X}) \propto p(\mathbf{y} \mid \mathbf{X}, \boldsymbol{\beta}, \sigma) \, p(\boldsymbol{\beta} \mid \sigma) \, p(\sigma) \] es decir: \[ p(\boldsymbol{\beta}, \sigma \mid \mathbf{y}, \mathbf{X}) \propto (\sigma^2)^{-n/2} \exp \left( -\frac{1}{2\sigma^2} (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^t(\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) \right) \exp \left( -\frac{1}{2}(\boldsymbol{\beta} - \boldsymbol{\mu}_0)^t \Sigma_0^{-1} (\boldsymbol{\beta}- \boldsymbol{\mu}_0) \right) (\sigma^2)^{-(a_0 + 1)} \exp \left( -\frac{b_0}{\sigma^2} \right) \]

Las distribuciones condicionales posterior de \(\boldsymbol{\beta}\) y \(\sigma^2\), se muestran a continuación.

Si la varianza \(\sigma^2\) es conocida (o condicionada a ella), la distribución posterior de \(\boldsymbol{\beta}\) está dada por

\[\boldsymbol{\beta} \mid \sigma, \mathbf{y}, \mathbf{X} \sim \mathcal{N}(\boldsymbol{\mu}_n, \boldsymbol{\Sigma}_n)\] Donde:

  • Media posterior: \[\boldsymbol{\mu}_n = \boldsymbol{\Sigma}_n\left(\boldsymbol{\Sigma}_0^{-1} \boldsymbol{\mu}_0+\frac{1}{\sigma^2}\mathbf{X}^t \mathbf{y} \right)\]

  • Matriz de precisión posterior: \[\boldsymbol{\Sigma}_n = \left(\boldsymbol{\Sigma}_0^{-1}+\frac{1}{\sigma^2}\mathbf{X}^t \mathbf{X}\right)^{-1}\] Por lo tanto: \[ p(\boldsymbol{\beta} \mid \sigma, \mathbf{y}, \mathbf{X}) = (2\pi)^{-\frac{k}{2}} |\boldsymbol{\Sigma}_n|^{-\frac{1}{2}} \exp\left( -\frac{1}{2} (\boldsymbol{\beta} - \boldsymbol{\mu}_n)^T \boldsymbol{\Sigma}_n^{-1} (\boldsymbol{\beta} - \boldsymbol{\mu}_n) \right) \]

Note que \[ p(\beta_0, \beta_1 \mid \sigma, \mathbf{y}, \mathbf{X}) = \frac{1}{2\pi \sigma_0 \sigma_1 \sqrt{1 - \rho^2}} \exp\left( -\frac{1}{2(1 - \rho^2)} \left[ \frac{(\beta_0 - \mu_0)^2}{\sigma_0^2} - \frac{2\rho(\beta_0 - \mu_0)(\beta_1 - \mu_1)}{\sigma_0 \sigma_1} + \frac{(\beta_1 - \mu_1)^2}{\sigma_1^2} \right] \right) \] de donde, se tiene: \[ p(\beta_0 \mid \sigma, \mathbf{y}, \mathbf{X}) = \frac{1}{\sqrt{2\pi \sigma_0^2}} \exp\left( -\frac{(\beta_0 - \mu_0)^2}{2\sigma_0^2} \right) \] \[ p(\beta_1 \mid \sigma, \mathbf{y}, \mathbf{X}) = \frac{1}{\sqrt{2\pi \sigma_1^2}} \exp\left( -\frac{(\beta_1 - \mu_1)^2}{2\sigma_1^2} \right) \] Distribución posterior marginal de la varianza \(\sigma\) está dada por:

\[\sigma \mid \boldsymbol{\beta}, \mathbf{y}, \mathbf{X} \sim \text{Inv-Gamma}(a_n, b_n)\] Donde

  • \(a_n = a_0 + \frac{n}{2}\)

  • \(b_n = b_0 + \frac{1}{2} \left( \mathbf{y} - X\boldsymbol{\beta})^t(\mathbf{y} - X\boldsymbol{\beta}\right)\)

\(n\) es el número de observaciones Por lo tanto: \[ p(\sigma \mid \boldsymbol{\beta}, \mathbf{y}, \mathbf{X}) = \frac{b_n^{a_n}}{\Gamma(a_n)} (\sigma^2)^{-(a_n + 1)} \exp\left( -\frac{b_n}{\sigma^2} \right) \]

Criterios de diagnóstico

Estas métricas describen la distribución de probabilidad posterior del parámetro:

Mean (Media Posterior): Es el primer momento de la distribución posterior. Representa el centro de masa de la probabilidad. Se interpreta como la estimación central (el valor más esperado).

  • Teórica (Continua): \[\mu = {E}[\theta | y] = \int_{-\infty}^{\infty} \theta \cdot p(\theta | y) \, d\theta\]

  • Muestral (MCMC): Si tienes \(N\) muestras guardadas (\(\theta_1, \theta_2, \dots, \theta_N\)), se calcula como el promedio aritmético simple: \[\bar{\theta} = \frac{1}{N}\sum_{i=1}^{N}\theta_i\]

StdDev (Desviación Estándar Posterior): Mide la dispersión cuadrática promedio alrededor de la media posterior. Una desviación estándar pequeña significa que la estimación es más precisa.

  • Teórica (Continua): \[\sigma=\sqrt{\text{Var}(\theta|y)}=\sqrt{\int_{-\infty}^{\infty}(\theta -\mu)^2 p(\theta|y) \,d\theta}\]

  • Muestral (MCMC): En la práctica, se utiliza la fórmula de la desviación estándar muestral corregida: \[s = \sqrt{\frac{1}{N-1} \sum_{i=1}^{N} (\theta_i - \bar{\theta})^2}\]

MAD (Desviación Absoluta de la Mediana Posterior): Es una medida robusta de dispersión basada en distancias absolutas en lugar de distancias al cuadrado. Es una medida de variabilidad similar a la desviación estándar, pero más robusta porque no se deja afectar tanto por valores extremos (outliers).

  • Muestral (MCMC): Primero se encuentra la mediana de las muestras (\(\tilde{\theta}\)). Luego, se calcula el valor absoluto de la desviación de cada muestra respecto a esa mediana, y finalmente se obtiene la mediana de esos valores absolutos. Para que sea comparable con la desviación estándar en una distribución normal, se multiplica por un factor de escala constante (\(\approx 1.4826\)): \[\text{MAD} = 1.4826 \times \text{mediana} \left( \left| \theta_i - \tilde{\theta} \right| \right)\] Donde \[\tilde{\theta} = \text{mediana}(\theta_1, \dots, \theta_N)\]

Mediana (Cuantil 50%): Es el valor que divide la masa de probabilidad exactamente en dos partes iguales, es decir, es el valor que parte las muestras a la mitad.

  • Teórica (Continua): Es el valor \(\tilde{\theta}\) tal que la función de distribución acumulada (CDF) es igual a 0.5: \[\int_{-\infty}^{\tilde{\theta}} p(\theta | y) \, d\theta = 0.5\]

  • Muestral (MCMC): Se ordenan las \(N\) muestras de menor a mayor (\(\theta_{(1)} \le \theta_{(2)} \le \dots \le \theta_{(N)}\)). Si \(N\) es impar, es el dato central; si \(N\) es par, es el promedio de los dos datos centrales: \[\tilde{\theta} = P_{50}\]

Intervalo de Credibilidad del 90% (Cuantiles 5% y 95%): A diferencia de la estadística frecuentista, aquí los límites describen una probabilidad directa sobre el parámetro. Este tipo de intervalo basado en cuantiles fijos se conoce como Intervalo de Igual Colas (Equal-Tailed Interval). El 5% y 95% delimitan un intervalo de credibilidad del 90%. Significa que hay un 90% de probabilidad de que el verdadero valor del parámetro se encuentre entre el cuantil 5% y el 95%.

  • Teórica (Continua): Buscamos los límites \(\theta_{0.05}\) y \(\theta_{0.95}\) tales que: \[\int_{-\infty}^{\theta_{0.05}} p(\theta | y) \, d\theta = 0.05 \quad \text{y} \quad \int_{-\infty}^{\theta_{0.95}} p(\theta | y) \, d\theta = 0.95\]

    Lo que garantiza que la masa central sea del 90 %: \[\int_{\theta_{0.05}}^{\theta_{0.95}} p(\theta | y) \, d\theta = 0.95 - 0.05 = 0.90\] Es decir, para el modelo que se tiene

  • Muestral (MCMC): Con las muestras ordenadas de menor a mayor, los límites corresponden a las posiciones de los cuantiles empíricos:

    • Límite inferior (\(\theta_{0.05}\)): Es el valor en la posición ordenada \(\lfloor 0.05 \times N \rfloor\).
    • Límite superior (\(\theta_{0.95}\)): Es el valor en la posición ordenada \(\lfloor 0.95 \times N \rfloor\).

Diagnósticos de Calidad del Muestreo:

R_hat (\(\hat{R}\)): Es el indicador de convergencia. Compara las diferentes cadenas de simulación para ver si todas llegaron al mismo resultado. Debe ser menor a 1.1. Si es mayor, significa que las cadenas no han convergido y los resultados no son válidos. 

El \(\hat{R}\) moderno no se calcula directamente sobre los valores crudos, sino transformando las muestras a rangos normales para corregir problemas con distribuciones de colas pesadas o asimétricas.

Para \(M\) cadenas y \(N\) simulaciones por cadena:

  1. Se transforman los valores originales \(G_{m,n}\) a rangos estandarizados \(U_{m,n}\) y luego a cuantiles normales: \[Z_{m,n} = \Phi^{-1}\left(\frac{U_{m,n} - 3/8}{M \cdot N + 1/4}\right)\] De esta forma se tendrá, las cadenas: Cadena 1: \[ Z_{1,1},...,Z_{1,N} \] Cadena 2: \[ Z_{2,1},...,Z_{2,N} \] Cadena M: \[ Z_{M,1},...,Z_{M,N} \]

  2. Con los cuantiles: Se calcula la varianza entre cadenas (\(B\)) y la varianza dentro de las cadenas, es decir, Donde:

  • \(W\) (Varianza dentro de cadenas): \[W = \frac{1}{M(N-1)} \sum_{m=1}^{M} \sum_{n=1}^{N} (Z_{m,n} - \bar{Z}_{m\cdot})^2\]
  • \(B\) (Varianza entre cadenas): \[B = \frac{N}{M-1} \sum_{m=1}^{M} (\bar{Z}_{m\cdot} - \bar{Z}_{\cdot\cdot})^2\] Dado como resultado el \(\hat{R}\): \[\hat{R} = \sqrt{\frac{W + \frac{1}{N}\left(B - W\right)}{W}}\]

ESS_bulk (Tamaño de Muestra Efectivo General): Mide la cantidad de muestras independientes equivalentes basándose en la autocorrelación de las cadenas en diferentes retrasos (lags), utilizando también la transformación de rangos normales (\(Z\)).

\[ESS_{bulk} = \frac{M \cdot N}{1 + 2 \sum_{k=1}^{\infty} \hat{\rho}_k}\]

Donde \(\hat{\rho}_k\) es la autocorrelación estimada al retraso \(k\) combinando todas las cadenas. La suma se detiene automáticamente en el retraso donde la autocorrelación deja de ser estadísticamente significativa (se vuelve negativa o ruidosa).

Por lo tanto, el ESS_bulk mide cuántas muestras verdaderamente independientes equivalen tus muestras de MCMC para estimar la media y la mediana. Debido a la autocorrelación, suele ser menor que el número total de pasos dados. Se recomienda que sea mayor a 400 por cadena. 

ESS_tail (Tamaño de Muestra Efectivo en las Colas): Para evaluar la precisión en los extremos de la distribución, se transforma la variable original a variables binarias (\(0\) o \(1\)) dependiendo de si caen o no en los cuantiles extremos (típicamente los percentiles del 5% y 95%, denotados como \(q_{\alpha}\)).

Se calcula un \(ESS\) individual para cada cuantil usando funciones indicadoras: \[I_{m,n} = \mathbb{I}(G_{m,n} \le q_{\alpha})\]

Luego se aplica la fórmula del ESS estándar sobre estos indicadores \(I_{m,n}\): \[ESS_{\alpha} = \frac{M \cdot N}{1 + 2 \sum_{k=1}^{\infty} \hat{\rho}_{k,\alpha}}\]

Finalmente, \(ESS_{tail}\) se define como el mínimo entre los tamaños de muestra efectivos de ambos extremos estudiados: \[ESS_{tail} = \min(ESS_{0.05},\, ESS_{0.95})\]

Así el ESS_tail mide la calidad de las muestras en las colas de la distribución (para calcular bien los percentiles como el 5% y 95%). También se recomienda que sea mayor a 400 por cadena. 

El MCSE (Error Estándar de Monte Carlo para la Media): Estima la desviación estándar teórica de tu estimación de la media debido al número finito de simulaciones y la dependencia entre ellas. Utiliza la desviación estándar de las muestras combinadas (\(\hat{\sigma}\)) dividida por la raíz cuadrada del tamaño de muestra efectivo:

\[MCSE(\mu) = \frac{\hat{\sigma}}{\sqrt{ESS_{bulk}}}\]

Donde \(\hat{\sigma}\) es la desviación estándar corregida de la población estimada a partir de las \(M \cdot N\) muestras.

Así el MCSE indica cuánta incertidumbre hay en la promedio debido a que usamos simulaciones en lugar de un cálculo exacto. Lo ideal es que el MCSE sea muy pequeño en comparación con la StdDev (regla general: menor al 5 % de la desviación estándar). 

ESS_bulk/s (Velocidad de Eficiencia Computacional): Es una métrica puramente aritmética de rendimiento de software/hardware. Divide el volumen de muestras independientes generadas para la media entre el tiempo de ejecución total en segundos exclusivamente de la fase de muestreo (\(t_{sampling}\)), excluyendo el warmup:

\[\text{ESS_bulk/s} = \frac{ESS_{bulk}}{t_{sampling}}\]

ESS_bulk/s: que mide la velocidad de eficiencia del algoritmo, indica cuántas muestras efectivas (“Bulk ESS”) está generando el modelo por cada segundo de reloj. Es útil para comparar el rendimiento computacional entre diferentes modelos. 

Modelo 4.4

ruta_dat<-"C:\\Documetos MGM_OS_Dell\\Documentos Maria_Dell\\Facultad_Matematicas\\1. Programa_Licenciatura en Matematicas\\Cursos de la LM\\Modelacion_bayesiana_LM\\Semestre_2_2025\\Bayesian_Statistical_Modeling_with_Stan_R_and_Python-master\\chap04\\input\\data-salary.csv" 

d <-read.csv(file=ruta_dat) # Datos
d <-read.csv(file=ruta_dat) # Datos


data<-list(N=nrow(d),X=d$X,Y=d$Y) 

model<- cmdstan_model(stan_file="C:\\Documetos MGM_OS_Dell\\Documentos Maria_Dell\\Facultad_Matematicas\\1. Programa_Licenciatura en Matematicas\\Cursos de la LM\\Modelacion_bayesiana_LM\\Semestre_2_2025\\Bayesian_Statistical_Modeling_with_Stan_R_and_Python-master\\chap04\\model\\model4-4.stan")

fit<-model$sample(data=data,
                  save_warmup = TRUE,
                  refresh         = 0        #  APAGA EL PROGRESO EN PANTALLA
                  )
## Running MCMC with 4 sequential chains...
## 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 1.2 seconds.

Resultados

a<-fit$summary()

a[2:4,] %>% 
  flextable() %>% 
  colformat_double(digits = 3) %>% 
  set_caption()

variable

mean

median

sd

mad

q5

q95

rhat

ess_bulk

ess_tail

a

38.738

38.712

1.481

1.418

36.370

41.161

1.004

1,303.983

1,311.453

b

0.751

0.752

0.091

0.088

0.600

0.903

1.004

1,222.605

1,578.463

sigma

2.905

2.810

0.630

0.578

2.071

4.086

1.001

1,610.661

1,676.688

pdf(file="C:\\Documetos MGM_OS_Dell\\Documentos Maria_Dell\\Facultad_Matematicas\\1. Programa_Licenciatura en Matematicas\\Cursos de la LM\\Modelacion_bayesiana_LM\\Semestre_2_2026\\cadenas.pdf")
plot(as_mcmc.list(fit)[,2:5]) 
dev.off()
## png 
##   2
# Get posterior draws
muestras<- fit$draws()  # Muestras de las distribuciones

as_draws_df(muestras)[,2:4] # Muestras de los parámetros

Histogramas de las cadenas

library(bayesplot) 
# 2. Guarda cada gráfico en una variable
p1 <- mcmc_hist(fit$draws("a"))
p2 <- mcmc_hist(fit$draws("b"))

p3 <- mcmc_hist(fit$draws("sigma"))

# En una sola columna (uno arriba del otro):
p1 / p2 / p3

Gráficos de las cadenas

 #inc_warmup = TRUE te permite ver las  cadenas desde sus puntos de inicio aleatorios (init_fun)
muest_par<- fit$draws(variables = c("a", "b", "sigma"))


t1 <- mcmc_trace(muest_par, pars = "a") + 
        ggtitle("Cadena de 'a'")
t2 <- mcmc_trace(muest_par, pars = "b") +
        ggtitle("Cadena de 'b'")
t3 <- mcmc_trace(muest_par, pars = "sigma") + 
        ggtitle("Cadena de 'sigma'")

# 3. Organizar los 3 gráficos en un solo arreglo (uno abajo del otro)
t1 / t2 / t3

Métricas del parámetro estimado

Para \[ \boldsymbol{\theta}=(\beta_0,\beta_1,\sigma^2) \] \(N\) muestras generadas: \(\boldsymbol{\theta}_1, \boldsymbol{\theta}_2, \dots, \boldsymbol{\theta}_N\), se tienen las muestras

(\(\beta_{01}, \beta_{02}, \dots, \beta_{0N}\)) (\(\beta_{11}, \beta_{12}, \dots, \beta_{1N}\)) (\(\sigma_1, \sigma_2, \dots, \sigma_N\)) Con las cuales se calculan las siguientes métricas.

  • Mean (Media Posterior): Estimación central de la distribución posterior

Si tiene \[\bar{\beta}_0 = \frac{1}{N}\sum_{i=1}^{N}\beta_{0i}\] \[\bar{\beta}_1 = \frac{1}{N}\sum_{i=1}^{N}\beta_{01}\] \[\bar{\sigma} = \frac{1}{N}\sum_{i=1}^{N}\sigma_i\] - StdDev (Desviación Estándar Posterior): Mide la dispersión cuadrática promedio alrededor de la media posterior. Una desviación estándar pequeña significa que la estimación es más precisa.

Entonces \[s_{\beta_0} = \sqrt{\frac{1}{N-1} \sum_{i=1}^{N} (\beta_{0i} - \bar{\beta}_0)^2}\] \[s_{\beta_1} = \sqrt{\frac{1}{N-1} \sum_{i=1}^{N} (\beta_{1i} - \bar{\beta}_1)^2}\] \[s_{\sigma } = \sqrt{\frac{1}{N-1} \sum_{i=1}^{N} (\sigma_{i} - \bar{\sigma})^2}\]

  • MAD (Desviación Absoluta de la Mediana Posterior):

Donde \[ \text{MAD} = 1.4826 \times \text{mediana} \left( \left| \theta_i - \tilde{\theta} \right| \right) \] Donde \[\tilde{\theta} = \text{mediana}(\theta_1, \dots, \theta_N)\]

  • Mediana (Cuantil 50%): Es el valor que divide la masa de probabilidad exactamente en dos partes iguales.

Es dicir, - Muestral (MCMC): Se ordenan las \(N\) muestras de menor a mayor (\(\theta_{(1)} \le \theta_{(2)} \le \dots \le \theta_{(N)}\)). Si \(N\) es impar, es el dato central; si \(N\) es par, es el promedio de los dos datos centrales: \[\tilde{\theta} = P_{50}\]

Para cada valor del -Intervalo de Credibilidad del 90% (Cuantiles 5% y 95%): El 5% y 95% delimitan un intervalo de credibilidad del 90%. Significa que hay un 90% de probabilidad de que el verdadero valor del parámetro se encuentre entre el cuantil 5% y el 95%.

Intervalo creible del 90 %: Con las muestras ordenadas de menor a mayor, los límites corresponden a las posiciones de los cuantiles empíricos: Límite inferior (\(\theta_{0.05}\)): Es el valor en la posición ordenada \[\lfloor 0.05 \times N \rfloor\]. Límite superior (\(\theta_{0.95}\)): Es el valor en la posición ordenada \[\lfloor 0.95 \times N \rfloor\].

Criterios de diagnóstico de las cadenas

  • R_hat (\(\hat{R}\)): Es el indicador de convergencia. Compara las diferentes cadenas de simulación para ver si todas llegaron al mismo resultado. Debe ser menor a 1.1. Si es mayor, significa que las cadenas no han convergido y los resultados no son válidos. 

  • ESS_bulk (Tamaño de Muestra Efectivo General): Mide cuántas muestras verdaderamente independientes equivalen tus muestras de MCMC para estimar la media y la mediana. Debido a la autocorrelación, suele ser menor que el número total de pasos dados. Se recomienda que sea mayor a 400 por cadena. 

  • ESS_tail (Tamaño de Muestra Efectivo en las Colas): Mide la calidad de las muestras en las colas de la distribución. También se recomienda que sea mayor a 400 por cadena. 

  • El MCSE (Error Estándar de Monte Carlo para la Media): Indica cuánta incertidumbre hay en la promedio debido a que usamos simulaciones en lugar de un cálculo exacto. Lo ideal es que el MCSE sea muy pequeño en comparación con la StdDev (regla general: menor al 5 % de la desviación estándar). 

  • ESS_bulk/s (Velocidad de Eficiencia Computacional): Es útil para comparar el rendimiento computacional entre diferentes modelos. 

\[\text{ESS_bulk/s} = \frac{ESS_{bulk}}{t_{sampling}}\] \(t_{sampling}\) excluye el warmup:

Modelo: 4.2.4 Adjust the Settings of MCMC

# Condiciones iniciales
init_fun<-function(chain_id)
          {  set.seed(chain_id) 
             list(a=runif(1,100,200),
                  b=runif(1,10,20),
                  sigma=rexp(1, rate = 20)  # valor positivo
                  ) 
            } 

fit<-model$sample(data=data,
                  seed=123,        # semilla
                  init=init_fun,   # condiciones iniciales
                  chains=2,        # Número de cadenas
                  iter_warmup=500,
                  iter_sampling=500,
                  thin=10,           # Adelgazamiento
                  parallel_chains=3, 
                  save_warmup=TRUE,  # para guardar los resultados de las iteraciones de calentamiento
                  refresh         = 0        #  APAGA EL PROGRESO EN PANTALLA
                  )
## Running MCMC with 2 chains, at most 3 in parallel...
## 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## 
## Both chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.

Muestras de los parámetros

d_ms<-fit$draws(format="df") # Muestras MCMC
quantile(d_ms$b, probs=c(0.025, 0.975)) # Intervalo creible de parámetro b
##      2.5%     97.5% 
## 0.5805591 0.9063826

Cada columna representa la destribución marginal del parámetro

\[ p(a|\mathbf{y}, \mathbf{X}) \] \[ p(b|\mathbf{y}, \mathbf{X}) \] \[ p(\sigma|\mathbf{y}, \mathbf{X}) \]

cadenas<-d_ms[1:10,2:4]

flextable(cadenas) %>% 
  set_caption("Primeros 10 valores de las cadenas de los parámetros")
Primeros 10 valores de las cadenas de los parámetros

a

b

sigma

37.91750

0.7955735

2.873143

38.78835

0.7499173

2.298703

37.45972

0.8103914

2.249689

38.00737

0.6923974

2.554942

37.97603

0.8485471

3.058927

37.26813

0.8549864

2.040516

39.77937

0.7200887

4.282367

37.02942

0.8010600

2.083628

39.49976

0.7157623

2.202287

38.30838

0.7621003

2.732307