Práctica de Análisis Bayesiano Multivariado aplicado a datos geofisicos del volcan Azufral
Gilbert Rodriguez y Diego Robayo
2026-05-16
En esta primera fase del análisis se realizó la carga y preparación de los datos para garantizar que cumplieran las condiciones necesarias del modelo bayesiano. Inicialmente, se importó el archivo Excel que contiene los 90 registros de campo utilizados en el estudio. Posteriormente, se aplicó un diagnóstico de normalidad para cada variable mediante una función que calculó estadísticos descriptivos básicos —como media, mediana, desviación estándar, valores mínimos y máximos— además de medidas de asimetría y curtosis, las cuales permiten evaluar qué tan cercana es la distribución de los datos a una distribución Normal. Complementariamente, se utilizaron las pruebas formales de Shapiro-Wilk y Lilliefors para contrastar la hipótesis nula de normalidad. Cuando alguna variable no cumplía este supuesto y todos sus valores eran positivos, se aplicó automáticamente una transformación logarítmica en base 10 con el fin de reducir sesgos y aproximar mejor la distribución a la Normal. Finalmente, todas las variables seleccionadas fueron estandarizadas utilizando la transformación z=(x−μ)/σ, logrando que cada una tuviera media cero y desviación estándar uno. Este paso fue fundamental debido a que las variables originales se encontraban en escalas físicas diferentes, como temperatura, resistividad y profundidad, y la estandarización permitió hacerlas comparables antes de incorporarlas al modelo bayesiano.
# -----------------------------------------------------------------------------
# 1. CARGA DE DATOS
# -----------------------------------------------------------------------------
# Ajusta la ruta si el archivo está en otra carpeta
df <- read_excel("D:/15.UNAL_Estadistica/TeoriaEstadisticaMultivariada/AnalisisBayesiano/SST_SEV_AZ.xlsx")
# Mostrar primeras filas y estructura
cat("ESTRUCTURA DEL DATASET")## ESTRUCTURA DEL DATASETprint(str(df))## tibble [90 × 8] (S3: tbl_df/tbl/data.frame)
## $ ID : num [1:90] 1 2 3 4 5 6 7 8 9 10 ...
## $ X_Magna : num [1:90] 935860 936818 938660 934697 935566 ...
## $ Y_Magna : num [1:90] 612035 611469 610491 611855 611458 ...
## $ Res_20cm : num [1:90] 148.3 16.3 85 251.9 440.9 ...
## $ T_20cm : num [1:90] 12.1 15.4 15.7 12.4 12.1 14.1 15.3 13 14.4 18 ...
## $ Res_150cm: num [1:90] 512 432 129 469 214 ...
## $ T_150cm : num [1:90] 12.5 15.7 16.4 12.1 12.5 14.8 15.5 12.9 14.1 17 ...
## $ Altura : num [1:90] 3212 3179 3006 3361 3212 ...
## NULLcat("PRIMERAS FILAS")## PRIMERAS FILASprint(head(df))## # A tibble: 6 × 8
## ID X_Magna Y_Magna Res_20cm T_20cm Res_150cm T_150cm Altura
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 935860 612035 148. 12.1 512. 12.5 3212.
## 2 2 936818 611469 16.3 15.4 432. 15.7 3179.
## 3 3 938660 610491 85.0 15.7 129. 16.4 3006
## 4 4 934697 611855 252. 12.4 469. 12.1 3361.
## 5 5 935566 611458 441. 12.1 214. 12.5 3212.
## 6 6 936073 610990 154. 14.1 112. 14.8 3086.cat("RESUMEN ESTADÍSTICO INICIAL")## RESUMEN ESTADÍSTICO INICIALprint(summary(df))## ID X_Magna Y_Magna Res_20cm
## Min. : 1.00 Min. :924352 Min. :599487 Min. : 15.72
## 1st Qu.: 24.25 1st Qu.:928846 1st Qu.:603789 1st Qu.: 110.26
## Median : 50.50 Median :933658 Median :605485 Median : 197.15
## Mean : 50.47 Mean :932727 Mean :605748 Mean : 291.86
## 3rd Qu.: 77.50 3rd Qu.:935947 3rd Qu.:607879 3rd Qu.: 339.90
## Max. :100.00 Max. :939391 Max. :612035 Max. :1744.40
## T_20cm Res_150cm T_150cm Altura
## Min. :12.10 Min. : 12.40 Min. :11.90 Min. :2909
## 1st Qu.:13.80 1st Qu.: 69.54 1st Qu.:14.03 1st Qu.:3013
## Median :14.30 Median : 163.02 Median :14.80 Median :3047
## Mean :14.92 Mean : 353.08 Mean :14.86 Mean :3056
## 3rd Qu.:15.60 3rd Qu.: 325.07 3rd Qu.:15.50 3rd Qu.:3086
## Max. :19.80 Max. :2390.00 Max. :18.90 Max. :3361# Seleccionar únicamente las 5 variables geofísicas de interés
# IMPORTANTE: ajusta los nombres de columna según tu archivo real
variables_originales <- c("T_20cm", "T_150cm", "Res_20cm", "Res_150cm", "Altura")
# Verificar que las columnas existan
cat("COLUMNAS DISPONIBLES EN EL ARCHIVO")## COLUMNAS DISPONIBLES EN EL ARCHIVOprint(names(df))## [1] "ID" "X_Magna" "Y_Magna" "Res_20cm" "T_20cm" "Res_150cm"
## [7] "T_150cm" "Altura"# Extraer subconjunto con las 5 variables
datos <- df %>%
dplyr::select(all_of(variables_originales))
cat("DATASET CON LAS 5 VARIABLES GEOFÍSICAS")## DATASET CON LAS 5 VARIABLES GEOFÍSICASprint(summary(datos))## T_20cm T_150cm Res_20cm Res_150cm
## Min. :12.10 Min. :11.90 Min. : 15.72 Min. : 12.40
## 1st Qu.:13.80 1st Qu.:14.03 1st Qu.: 110.26 1st Qu.: 69.54
## Median :14.30 Median :14.80 Median : 197.15 Median : 163.02
## Mean :14.92 Mean :14.86 Mean : 291.86 Mean : 353.08
## 3rd Qu.:15.60 3rd Qu.:15.50 3rd Qu.: 339.90 3rd Qu.: 325.07
## Max. :19.80 Max. :18.90 Max. :1744.40 Max. :2390.00
## Altura
## Min. :2909
## 1st Qu.:3013
## Median :3047
## Mean :3056
## 3rd Qu.:3086
## Max. :3361# -----------------------------------------------------------------------------
# 2. FUNCIÓN AUXILIAR: DIAGNÓSTICO DE NORMALIDAD POR VARIABLE
# -----------------------------------------------------------------------------
diagnostico_normalidad <- function(x, nombre_var) {
cat("\n", rep("=", 55), "\n", sep = "")
cat(" VARIABLE:", nombre_var, "\n")
cat(rep("=", 55), "\n", sep = "")
# --- Estadísticos descriptivos ---
cat(" Media :", round(mean(x, na.rm = TRUE), 4), "\n")
cat(" Mediana :", round(median(x, na.rm = TRUE), 4), "\n")
cat(" Desv. estándar :", round(sd(x, na.rm = TRUE), 4), "\n")
cat(" Mínimo :", round(min(x, na.rm = TRUE), 4), "\n")
cat(" Máximo :", round(max(x, na.rm = TRUE), 4), "\n")
# --- Asimetría y curtosis ---
asim <- skewness(x, na.rm = TRUE)
curt <- kurtosis(x, na.rm = TRUE) - 3 # curtosis excesiva (0 en Normal)
cat(" Asimetría :", round(asim, 4),
ifelse(abs(asim) < 0.5, " -> Aproximadamente simétrica",
ifelse(abs(asim) < 1, " -> Asimetría moderada",
" -> Asimetría fuerte")), "\n")
cat(" Curtosis exc. :", round(curt, 4),
ifelse(abs(curt) < 1, " -> Curtosis aceptable",
" -> Curtosis elevada (colas pesadas)"), "\n")
# --- Pruebas formales de normalidad ---
# Shapiro-Wilk (mejor para n < 50; aún válida con n <= 5000)
sw <- shapiro.test(x)
# Kolmogorov-Smirnov con corrección de Lilliefors
lil <- lillie.test(x)
cat("\n >> Prueba Shapiro-Wilk:\n")
cat(" W =", round(sw$statistic, 4),
" p-value =", round(sw$p.value, 5), "\n")
cat(" Conclusión:",
ifelse(sw$p.value > 0.05,
"No se rechaza H0 -> distribución compatible con Normal",
"Se rechaza H0 -> distribución NO es Normal (p < 0.05)"), "\n")
cat("\n >> Prueba Lilliefors (KS corregido):\n")
cat(" D =", round(lil$statistic, 4),
" p-value =", round(lil$p.value, 5), "\n")
cat(" Conclusión:",
ifelse(lil$p.value > 0.05,
"No se rechaza H0 -> distribución compatible con Normal",
"Se rechaza H0 -> distribución NO es Normal (p < 0.05)"), "\n")
# --- Decisión de transformación ---
necesita_transformacion <- (sw$p.value < 0.05 | lil$p.value < 0.05) &
all(x > 0, na.rm = TRUE)
cat("\n >> DECISIÓN DE TRANSFORMACIÓN:\n")
if (necesita_transformacion) {
cat(" ⚠ La variable NO sigue distribución Normal.\n")
cat(" Todos los valores son positivos -> se aplicará log10(x).\n")
# Verificar si la transformación ayuda
x_log <- log10(x)
sw_log <- shapiro.test(x_log)
cat(" Shapiro-Wilk sobre log10(", nombre_var, "): p =",
round(sw_log$p.value, 5), "\n")
if (sw_log$p.value > 0.05) {
cat(" ✓ Transformación log10 NORMALIZA la variable.\n")
} else {
cat(" ✗ Transformación log10 mejora pero no normaliza completamente.\n")
cat(" Se usará de todas formas (práctica estándar en geofísica).\n")
}
} else if (sw$p.value < 0.05 | lil$p.value < 0.05) {
cat(" ⚠ La variable NO sigue distribución Normal.\n")
cat(" Tiene valores <= 0 -> NO se puede aplicar log10.\n")
cat(" Se usará la variable original (verificar outliers).\n")
} else {
cat(" ✓ La variable es compatible con distribución Normal.\n")
cat(" No se requiere transformación.\n")
}
return(list(
necesita_log = necesita_transformacion,
p_sw = sw$p.value,
p_lil = lil$p.value,
asimetria = asim
))
}
# -----------------------------------------------------------------------------
# 3. DIAGNÓSTICO DE NORMALIDAD PARA LAS 5 VARIABLES
# -----------------------------------------------------------------------------
cat("DIAGNÓSTICO DE NORMALIDAD - VARIABLES ORIGINALES")## DIAGNÓSTICO DE NORMALIDAD - VARIABLES ORIGINALESresultados_diag <- lapply(variables_originales, function(v) {
diagnostico_normalidad(datos[[v]], v)
})##
## =======================================================
## VARIABLE: T_20cm
## =======================================================
## Media : 14.9189
## Mediana : 14.3
## Desv. estándar : 1.6395
## Mínimo : 12.1
## Máximo : 19.8
## Asimetría : 0.9802 -> Asimetría moderada
## Curtosis exc. : 0.5506 -> Curtosis aceptable
##
## >> Prueba Shapiro-Wilk:
## W = 0.9064 p-value = 1e-05
## Conclusión: Se rechaza H0 -> distribución NO es Normal (p < 0.05)
##
## >> Prueba Lilliefors (KS corregido):
## D = 0.1642 p-value = 0
## Conclusión: Se rechaza H0 -> distribución NO es Normal (p < 0.05)
##
## >> DECISIÓN DE TRANSFORMACIÓN:
## ⚠ La variable NO sigue distribución Normal.
## Todos los valores son positivos -> se aplicará log10(x).
## Shapiro-Wilk sobre log10( T_20cm ): p = 0.00014
## ✗ Transformación log10 mejora pero no normaliza completamente.
## Se usará de todas formas (práctica estándar en geofísica).
##
## =======================================================
## VARIABLE: T_150cm
## =======================================================
## Media : 14.8622
## Mediana : 14.8
## Desv. estándar : 1.3763
## Mínimo : 11.9
## Máximo : 18.9
## Asimetría : 0.3804 -> Aproximadamente simétrica
## Curtosis exc. : 0.4713 -> Curtosis aceptable
##
## >> Prueba Shapiro-Wilk:
## W = 0.9771 p-value = 0.11223
## Conclusión: No se rechaza H0 -> distribución compatible con Normal
##
## >> Prueba Lilliefors (KS corregido):
## D = 0.0877 p-value = 0.08404
## Conclusión: No se rechaza H0 -> distribución compatible con Normal
##
## >> DECISIÓN DE TRANSFORMACIÓN:
## ✓ La variable es compatible con distribución Normal.
## No se requiere transformación.
##
## =======================================================
## VARIABLE: Res_20cm
## =======================================================
## Media : 291.8598
## Mediana : 197.15
## Desv. estándar : 310.5223
## Mínimo : 15.72
## Máximo : 1744.4
## Asimetría : 2.5913 -> Asimetría fuerte
## Curtosis exc. : 7.6273 -> Curtosis elevada (colas pesadas)
##
## >> Prueba Shapiro-Wilk:
## W = 0.7101 p-value = 0
## Conclusión: Se rechaza H0 -> distribución NO es Normal (p < 0.05)
##
## >> Prueba Lilliefors (KS corregido):
## D = 0.1959 p-value = 0
## Conclusión: Se rechaza H0 -> distribución NO es Normal (p < 0.05)
##
## >> DECISIÓN DE TRANSFORMACIÓN:
## ⚠ La variable NO sigue distribución Normal.
## Todos los valores son positivos -> se aplicará log10(x).
## Shapiro-Wilk sobre log10( Res_20cm ): p = 0.80959
## ✓ Transformación log10 NORMALIZA la variable.
##
## =======================================================
## VARIABLE: Res_150cm
## =======================================================
## Media : 353.0823
## Mediana : 163.0245
## Desv. estándar : 487.4336
## Mínimo : 12.4
## Máximo : 2390
## Asimetría : 2.3246 -> Asimetría fuerte
## Curtosis exc. : 5.1259 -> Curtosis elevada (colas pesadas)
##
## >> Prueba Shapiro-Wilk:
## W = 0.6643 p-value = 0
## Conclusión: Se rechaza H0 -> distribución NO es Normal (p < 0.05)
##
## >> Prueba Lilliefors (KS corregido):
## D = 0.278 p-value = 0
## Conclusión: Se rechaza H0 -> distribución NO es Normal (p < 0.05)
##
## >> DECISIÓN DE TRANSFORMACIÓN:
## ⚠ La variable NO sigue distribución Normal.
## Todos los valores son positivos -> se aplicará log10(x).
## Shapiro-Wilk sobre log10( Res_150cm ): p = 0.33398
## ✓ Transformación log10 NORMALIZA la variable.
##
## =======================================================
## VARIABLE: Altura
## =======================================================
## Media : 3056.45
## Mediana : 3046.823
## Desv. estándar : 77.9917
## Mínimo : 2908.624
## Máximo : 3360.951
## Asimetría : 1.032 -> Asimetría fuerte
## Curtosis exc. : 1.8992 -> Curtosis elevada (colas pesadas)
##
## >> Prueba Shapiro-Wilk:
## W = 0.938 p-value = 0.00033
## Conclusión: Se rechaza H0 -> distribución NO es Normal (p < 0.05)
##
## >> Prueba Lilliefors (KS corregido):
## D = 0.1213 p-value = 0.00232
## Conclusión: Se rechaza H0 -> distribución NO es Normal (p < 0.05)
##
## >> DECISIÓN DE TRANSFORMACIÓN:
## ⚠ La variable NO sigue distribución Normal.
## Todos los valores son positivos -> se aplicará log10(x).
## Shapiro-Wilk sobre log10( Altura ): p = 0.00103
## ✗ Transformación log10 mejora pero no normaliza completamente.
## Se usará de todas formas (práctica estándar en geofísica).names(resultados_diag) <- variables_originales
# -----------------------------------------------------------------------------
# 4. APLICAR TRANSFORMACIONES DONDE SEA NECESARIO
# -----------------------------------------------------------------------------
cat("APLICANDO TRANSFORMACIONES")## APLICANDO TRANSFORMACIONESdatos_transf <- datos # copia de trabajo
for (v in variables_originales) {
if (resultados_diag[[v]]$necesita_log) {
nuevo_nombre <- paste0("log10_", v)
datos_transf[[nuevo_nombre]] <- log10(datos_transf[[v]])
cat(" -> Se creó la columna", nuevo_nombre,
"como log10(", v, ")\n")
}
}## -> Se creó la columna log10_T_20cm como log10( T_20cm )
## -> Se creó la columna log10_Res_20cm como log10( Res_20cm )
## -> Se creó la columna log10_Res_150cm como log10( Res_150cm )
## -> Se creó la columna log10_Altura como log10( Altura )# Identificar qué columnas usar para el análisis (originales o transformadas)
cols_analisis <- sapply(variables_originales, function(v) {
if (resultados_diag[[v]]$necesita_log) paste0("log10_", v) else v
})
cat("Variables que entrarán al análisis bayesiano:\n")## Variables que entrarán al análisis bayesiano:for (i in seq_along(variables_originales)) {
cat(" ", variables_originales[i], "->", cols_analisis[i], "\n")
}## T_20cm -> log10_T_20cm
## T_150cm -> T_150cm
## Res_20cm -> log10_Res_20cm
## Res_150cm -> log10_Res_150cm
## Altura -> log10_Altura# -----------------------------------------------------------------------------
# 5. ESTANDARIZACIÓN (media = 0, desviación estándar = 1)
# -----------------------------------------------------------------------------
cat("ESTANDARIZACIÓN DE VARIABLES")## ESTANDARIZACIÓN DE VARIABLEScat("Fórmula: z = (x - media) / desviación_estándar\n\n")## Fórmula: z = (x - media) / desviación_estándardatos_std <- as.data.frame(scale(datos_transf[, cols_analisis]))
# Renombrar columnas estandarizadas con prefijo "std_"
colnames(datos_std) <- paste0("std_", variables_originales)
cat("Medias después de estandarizar (deben ser ~0):\n")## Medias después de estandarizar (deben ser ~0):print(round(colMeans(datos_std), 6))## std_T_20cm std_T_150cm std_Res_20cm std_Res_150cm std_Altura
## 0 0 0 0 0cat("Desviaciones estándar después de estandarizar (deben ser ~1):\n")## Desviaciones estándar después de estandarizar (deben ser ~1):print(round(apply(datos_std, 2, sd), 6))## std_T_20cm std_T_150cm std_Res_20cm std_Res_150cm std_Altura
## 1 1 1 1 1cat("Resumen de las variables estandarizadas:\n")## Resumen de las variables estandarizadas:print(summary(datos_std))## std_T_20cm std_T_150cm std_Res_20cm std_Res_150cm
## Min. :-1.9300 Min. :-2.15226 Min. :-2.64569 Min. :-2.16513
## 1st Qu.:-0.6849 1st Qu.:-0.60830 1st Qu.:-0.58370 1st Qu.:-0.74363
## Median :-0.3478 Median :-0.04521 Median : 0.03196 Median :-0.04091
## Mean : 0.0000 Mean : 0.00000 Mean : 0.00000 Mean : 0.00000
## 3rd Qu.: 0.4764 3rd Qu.: 0.46339 3rd Qu.: 0.60902 3rd Qu.: 0.52816
## Max. : 2.7345 Max. : 2.93373 Max. : 2.34087 Max. : 2.17313
## std_Altura
## Min. :-1.9536
## 1st Qu.:-0.5569
## Median :-0.1125
## Mean : 0.0000
## 3rd Qu.: 0.3900
## Max. : 3.7791# -----------------------------------------------------------------------------
# FUNCIÓN AUXILIAR PARA HISTOGRAMAS
# -----------------------------------------------------------------------------
crear_histograma <- function(x, titulo,
color = "#2E86AB",
binwidth = NULL) {
df_plot <- data.frame(valor = x)
p <- ggplot(df_plot, aes(x = valor)) +
geom_histogram(
aes(y = after_stat(density)),
fill = color,
color = "white",
alpha = 0.8,
bins = ifelse(is.null(binwidth), 15, NULL),
binwidth = binwidth
) +
geom_density(
color = "#E63946",
linewidth = 1.0
) +
stat_function(
fun = dnorm,
args = list(
mean = mean(x, na.rm = TRUE),
sd = sd(x, na.rm = TRUE)
),
color = "#2D6A4F",
linewidth = 0.9,
linetype = "dashed"
) +
labs(
title = titulo,
subtitle = paste0(
"Asim: ",
round(skewness(x, na.rm = TRUE), 3),
" | Kurt: ",
round(kurtosis(x, na.rm = TRUE) - 3, 3)
),
x = "Valor",
y = "Densidad"
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(
face = "bold",
size = 12
),
plot.subtitle = element_text(
size = 10,
color = "gray40"
)
)
return(p)
}
# -----------------------------------------------------------------------------
# FIGURA 1A — TEMPERATURAS
# -----------------------------------------------------------------------------
plots_temp <- list(
crear_histograma(
datos[["T_20cm"]],
"T_20cm (original)",
color = "#2E86AB"
),
crear_histograma(
datos_transf[["log10_T_20cm"]],
"log10(T_20cm) — transformada",
color = "#F4A261"
),
crear_histograma(
datos[["T_150cm"]],
"T_150cm (original)",
color = "#2E86AB"
)
)
# Organización tipo matriz 2x2
fig1A <- grid.arrange(
grobs = plots_temp,
ncol = 2,
nrow = 2,
top = "FIGURA 1A — Variables de Temperatura"
)
# Guardar figura
ggsave(
"fig1A_temperaturas.png",
fig1A,
width = 12,
height = 8,
dpi = 250,
bg = "white"
)
# -----------------------------------------------------------------------------
# FIGURA 1B — RESISTIVIDADES
# -----------------------------------------------------------------------------
plots_res <- list(
crear_histograma(
datos[["Res_20cm"]],
"Res_20cm (original)",
color = "#2E86AB"
),
crear_histograma(
datos_transf[["log10_Res_20cm"]],
"log10(Res_20cm) — transformada",
color = "#F4A261"
),
crear_histograma(
datos[["Res_150cm"]],
"Res_150cm (original)",
color = "#2E86AB"
),
crear_histograma(
datos_transf[["log10_Res_150cm"]],
"log10(Res_150cm) — transformada",
color = "#F4A261"
)
)
# Organizar en matriz 2x2
fig1B <- grid.arrange(
grobs = plots_res,
ncol = 2,
nrow = 2,
top = "FIGURA 1B — Variables de Resistividad"
)
# Guardar figura
ggsave(
"fig1B_resistividades.png",
fig1B,
width = 12,
height = 8,
dpi = 250,
bg = "white"
)
# -----------------------------------------------------------------------------
# FIGURA 1C — ALTURA
# -----------------------------------------------------------------------------
plots_alt <- list(
crear_histograma(
datos[["Altura"]],
"Altura (original)",
color = "#2E86AB"
),
crear_histograma(
datos_transf[["log10_Altura"]],
"log10(Altura) — transformada",
color = "#F4A261"
)
)
fig1C <- grid.arrange(
grobs = plots_alt,
ncol = 1,
top = "FIGURA 1C — Altura"
)
ggsave(
"fig1C_altura.png",
fig1C,
width = 8,
height = 6,
dpi = 250,
bg = "white"
)
# -----------------------------------------------------------------------------
# 9. FIGURA 4: TABLA RESUMEN DE DIAGNÓSTICO
# -----------------------------------------------------------------------------
tabla_resumen <- data.frame(
Variable = variables_originales,
Media = round(sapply(variables_originales, function(v) mean(datos[[v]], na.rm = TRUE)), 3),
SD = round(sapply(variables_originales, function(v) sd(datos[[v]], na.rm = TRUE)), 3),
Asimetria = round(sapply(variables_originales, function(v) skewness(datos[[v]], na.rm = TRUE)), 3),
p_Shapiro = round(sapply(variables_originales, function(v) resultados_diag[[v]]$p_sw), 4),
p_Lilliefors = round(sapply(variables_originales, function(v) resultados_diag[[v]]$p_lil), 4),
Normal = sapply(variables_originales,
function(v) ifelse(resultados_diag[[v]]$p_sw > 0.05 &
resultados_diag[[v]]$p_lil > 0.05,
"Sí", "No")),
Transformada = sapply(variables_originales,
function(v) ifelse(resultados_diag[[v]]$necesita_log,
paste0("log10(", v, ")"), "Sin cambio")),
stringsAsFactors = FALSE
)
cat("TABLA RESUMEN DEL DIAGNÓSTICO DE NORMALIDAD")## TABLA RESUMEN DEL DIAGNÓSTICO DE NORMALIDADprint(tabla_resumen)## Variable Media SD Asimetria p_Shapiro p_Lilliefors Normal
## T_20cm T_20cm 14.919 1.639 0.980 0.0000 0.0000 No
## T_150cm T_150cm 14.862 1.376 0.380 0.1122 0.0840 Sí
## Res_20cm Res_20cm 291.860 310.522 2.591 0.0000 0.0000 No
## Res_150cm Res_150cm 353.082 487.434 2.325 0.0000 0.0000 No
## Altura Altura 3056.450 77.992 1.032 0.0003 0.0023 No
## Transformada
## T_20cm log10(T_20cm)
## T_150cm Sin cambio
## Res_20cm log10(Res_20cm)
## Res_150cm log10(Res_150cm)
## Altura log10(Altura)fig4 <- tableGrob(tabla_resumen, rows = NULL,
theme = ttheme_minimal(
core = list(fg_params = list(fontsize = 8)),
colhead = list(fg_params = list(fontsize = 8, fontface = "bold"),
bg_params = list(fill = "#2E86AB", col = NA, alpha = 0.8))
))
titulo_tabla <- grid::textGrob(
"FIGURA 4 — Tabla Resumen: Diagnóstico de Normalidad y Transformaciones Aplicadas",
gp = grid::gpar(fontface = "bold", fontsize = 10)
)
fig4_final <- gridExtra::arrangeGrob(titulo_tabla, fig4, ncol = 1,
heights = c(0.1, 0.9))
ggsave("fig4_tabla_diagnostico.png", fig4_final,
width = 12, height = 3, dpi = 180, bg = "white")
# -----------------------------------------------------------------------------
# 10. GUARDAR EL DATASET ESTANDARIZADO PARA LA FASE 2
# -----------------------------------------------------------------------------
# Exportar como CSV para usar en la siguiente fase del análisis bayesiano
write.csv(datos_std, "datos_estandarizados_fase2.csv", row.names = FALSE)En esta segunda fase se definió la estructura del modelo bayesiano y las distribuciones a priori de sus parámetros. Primero, se establecieron las dimensiones del problema mediante el número de observaciones (n=90) y el número de variables analizadas (p=5), valores que determinan la estructura matemática del modelo. Para el vector de medias μ, se especificó una distribución Normal multivariada con media a priori μ=(0,0,0,0,0) y matriz de covarianza Λ0=10I. Dado que las variables ya habían sido estandarizadas, asumir medias cercanas a cero resulta coherente, mientras que la matriz diagonal con valores altos representa una prior difusa o poco informativa, permitiendo que los datos tengan mayor influencia en la estimación posterior. Por otra parte, la matriz de covarianza Σ se modeló mediante una distribución Inverse-Wishart con grados de libertad v0=p+2=7 y matriz de escala S0=I. Esta distribución se utiliza de forma estándar en modelos bayesianos multivariados porque garantiza que la matriz de covarianza sea positiva definida. Además, la elección de pocos grados de libertad y de la matriz identidad como escala refleja nuevamente una prior débil, sin asumir previamente una estructura específica de correlación o variabilidad entre las variables.
# -----------------------------------------------------------------------------
# 1. CARGAR DATASET ESTANDARIZADO
# -----------------------------------------------------------------------------
datos_std <- read.csv("datos_estandarizados_fase2.csv")
cat("DATASET ESTANDARIZADO - FASE 2")## DATASET ESTANDARIZADO - FASE 2print(head(datos_std))## std_T_20cm std_T_150cm std_Res_20cm std_Res_150cm std_Altura
## 1 -1.9300360 -1.71631835 -0.2693755 0.9025389 1.9846967
## 2 0.3541730 0.60870368 -2.6084272 0.7635543 1.5722185
## 3 0.5369120 1.11730225 -0.8584069 -0.2344461 -0.6475368
## 4 -1.6980650 -2.00694610 0.2918347 0.8301781 3.7791024
## 5 -1.9300360 -1.71631835 0.8845800 0.1815980 1.9846967
## 6 -0.4811599 -0.04520876 -0.2300936 -0.3520920 0.3944537print(summary(datos_std))## std_T_20cm std_T_150cm std_Res_20cm std_Res_150cm
## Min. :-1.9300 Min. :-2.15226 Min. :-2.64569 Min. :-2.16513
## 1st Qu.:-0.6849 1st Qu.:-0.60830 1st Qu.:-0.58370 1st Qu.:-0.74363
## Median :-0.3478 Median :-0.04521 Median : 0.03196 Median :-0.04091
## Mean : 0.0000 Mean : 0.00000 Mean : 0.00000 Mean : 0.00000
## 3rd Qu.: 0.4764 3rd Qu.: 0.46339 3rd Qu.: 0.60902 3rd Qu.: 0.52816
## Max. : 2.7345 Max. : 2.93373 Max. : 2.34087 Max. : 2.17313
## std_Altura
## Min. :-1.9536
## 1st Qu.:-0.5569
## Median :-0.1125
## Mean : 0.0000
## 3rd Qu.: 0.3900
## Max. : 3.7791# -----------------------------------------------------------------------------
# 2. DIMENSIONES DEL PROBLEMA
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# 2. DIMENSIONES DEL PROBLEMA
# -----------------------------------------------------------------------------
n <- nrow(datos_std)
p <- ncol(datos_std)
cat("DIMENSIONES DEL PROBLEMA\n")## DIMENSIONES DEL PROBLEMAtabla_dim <- data.frame(
Concepto = c("Número de observaciones (n)",
"Número de variables (p)"),
Valor = c(n, p)
)
print(tabla_dim)## Concepto Valor
## 1 Número de observaciones (n) 90
## 2 Número de variables (p) 5# -----------------------------------------------------------------------------
# 2.1 PARÁMETROS CIENTÍFICOS DE INTERÉS
# -----------------------------------------------------------------------------
# Objetivo científico:
#
# 1) rho(T20, Res20)
# 2) rho(T150, Res150)
#
# Estas correlaciones se derivan de Sigma.
# -----------------------------------------------------------------------------
cat("PARÁMETROS CIENTÍFICOS DE INTERÉS\n")## PARÁMETROS CIENTÍFICOS DE INTERÉS# -----------------------------------------------------------------------------
# IDENTIFICAR ÍNDICES DE VARIABLES
# -----------------------------------------------------------------------------
idx_T20 <- which(colnames(datos_std) == "std_T_20cm")
idx_Res20 <- which(colnames(datos_std) == "std_Res_20cm")
idx_T150 <- which(colnames(datos_std) == "std_T_150cm")
idx_Res150 <- which(colnames(datos_std) == "std_Res_150cm")
idx_Alt <- which(colnames(datos_std) == "std_Altura")
# Tabla organizada de variables
tabla_variables <- data.frame(
Variable = c("Temperatura 20 cm",
"Resistividad 20 cm",
"Temperatura 150 cm",
"Resistividad 150 cm"),
Nombre_Columna = c("std_T_20cm",
"std_Res_20cm",
"std_T_150cm",
"std_Res_150cm"),
Indice = c(idx_T20,
idx_Res20,
idx_T150,
idx_Res150)
)
cat("Tabla de variables e índices:\n")## Tabla de variables e índices:print(tabla_variables)## Variable Nombre_Columna Indice
## 1 Temperatura 20 cm std_T_20cm 1
## 2 Resistividad 20 cm std_Res_20cm 3
## 3 Temperatura 150 cm std_T_150cm 2
## 4 Resistividad 150 cm std_Res_150cm 4# -----------------------------------------------------------------------------
# CORRELACIONES MUESTRALES INICIALES
# -----------------------------------------------------------------------------
# Estas NO son las correlaciones bayesianas finales.
# Son únicamente estimaciones exploratorias iniciales.
# -----------------------------------------------------------------------------
rho_init_T20_Res20 <- cor(
datos_std[, idx_T20],
datos_std[, idx_Res20]
)
rho_init_T150_Res150 <- cor(
datos_std[, idx_T150],
datos_std[, idx_Res150]
)
# Tabla organizada de correlaciones
tabla_correlaciones <- data.frame(
Parametro = c("rho(T20, Res20)",
"rho(T150, Res150)"),
Correlacion_Muestral = c(
rho_init_T20_Res20,
rho_init_T150_Res150
)
)
# Redondear SOLO columnas numéricas
tabla_correlaciones$Correlacion_Muestral <-
round(tabla_correlaciones$Correlacion_Muestral, 4)
cat("CORRELACIONES MUESTRALES INICIALES\n")## CORRELACIONES MUESTRALES INICIALESprint(tabla_correlaciones)## Parametro Correlacion_Muestral
## 1 rho(T20, Res20) -0.3136
## 2 rho(T150, Res150) -0.1974# -----------------------------------------------------------------------------
# 3. DEFINICIÓN DEL MODELO MULTIVARIANTE
# -----------------------------------------------------------------------------
# Supuesto:
#
# Y_i ~ MVN(mu, Sigma)
#
# donde:
#
# mu -> vector de medias (p x 1)
# Sigma -> matriz de covarianza (p x p)
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# 4. PRIORS PARA EL VECTOR DE MEDIAS mu
# -----------------------------------------------------------------------------
# Como los datos están estandarizados:
#
# media ≈ 0
# sd ≈ 1
#
# entonces usamos priors débiles centradas en 0.
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# 4. PRIOR PARA EL VECTOR DE MEDIAS mu
# -----------------------------------------------------------------------------
# Como los datos están estandarizados:
#
# media ≈ 0
# sd ≈ 1
#
# usamos una prior débil centrada en 0.
# -----------------------------------------------------------------------------
# ─── AQUÍ EMPIEZA LO NUEVO ───────────────────────────────────────────
# Índices de variables (ya los tienes definidos arriba, no los repitas)
# idx_T20, idx_T150, idx_Res20, idx_Res150, idx_Alt
# ya están definidos en tu código — no los vuelvas a definir
# PRIOR PARA mu
n0 <- 40
mu0 <- rep(0, p)
mu0[idx_T20] <- 0.30
mu0[idx_T150] <- 0.20
mu0[idx_Res20] <- -0.20
mu0[idx_Res150] <- -0.20
mu0[idx_Alt] <- 0.00
# Incertidumbre sobre mu0
Lambda0 <- diag(1 / n0, p)Distribucion prior para el vector de medias:
\[ \mu \sim N_p(\mu_0, \Lambda_0) \]
# -----------------------------------------------------------------------------
# TABLA DEL VECTOR mu0
# -----------------------------------------------------------------------------
tabla_mu0 <- data.frame(
Variable = colnames(datos_std),
Media_Prior = mu0
)
cat("Tabla del vector de medias a priori (mu0):\n")## Tabla del vector de medias a priori (mu0):print(tabla_mu0)## Variable Media_Prior
## 1 std_T_20cm 0.3
## 2 std_T_150cm 0.2
## 3 std_Res_20cm -0.2
## 4 std_Res_150cm -0.2
## 5 std_Altura 0.0# -----------------------------------------------------------------------------
# TABLA DE LA MATRIZ Lambda0
# -----------------------------------------------------------------------------
Lambda0_tabla <- as.data.frame(Lambda0)
colnames(Lambda0_tabla) <- colnames(datos_std)
rownames(Lambda0_tabla) <- colnames(datos_std)
cat("Matriz de covarianza prior Lambda0:\n")## Matriz de covarianza prior Lambda0:print(round(Lambda0_tabla, 4))## std_T_20cm std_T_150cm std_Res_20cm std_Res_150cm std_Altura
## std_T_20cm 0.025 0.000 0.000 0.000 0.000
## std_T_150cm 0.000 0.025 0.000 0.000 0.000
## std_Res_20cm 0.000 0.000 0.025 0.000 0.000
## std_Res_150cm 0.000 0.000 0.000 0.025 0.000
## std_Altura 0.000 0.000 0.000 0.000 0.025# -----------------------------------------------------------------------------
# 5. PRIOR PARA LA MATRIZ DE COVARIANZA Sigma
# -----------------------------------------------------------------------------
# Usaremos una distribución Inverse-Wishart:
#
# Sigma ~ IW(v0, S0)
#
# -----------------------------------------------------------------------------
# Grados de libertad
v0 <- p + 2
# Matriz de escala
S0 <- diag(1, p)Distribucion prior para la matriz de covarianza:
\[ \Sigma \sim \text{Inverse-Wishart}(v_0, S_0) \]
# chunk Fase2_3 — estructura que debes tener
# PRIOR PARA Sigma
S0 <- diag(1, p)
rownames(S0) <- colnames(datos_std)
colnames(S0) <- colnames(datos_std)
# Correlaciones a priori
S0[idx_Alt, idx_T20] <- -0.40; S0[idx_T20, idx_Alt] <- -0.40
S0[idx_Alt, idx_T150] <- -0.35; S0[idx_T150, idx_Alt] <- -0.35
S0[idx_T20, idx_Res20] <- -0.30; S0[idx_Res20, idx_T20] <- -0.30
S0[idx_T150, idx_Res150]<- -0.30; S0[idx_Res150,idx_T150] <- -0.30
S0[idx_T20, idx_T150] <- 0.50; S0[idx_T150, idx_T20] <- 0.50
S0[idx_Res20,idx_Res150]<- 0.40; S0[idx_Res150,idx_Res20]<- 0.40
# Verificar positiva definida — OBLIGATORIO
autovalores <- eigen(S0)$values
cat("Autovalores de S0:", round(autovalores, 4), "\n")## Autovalores de S0: 1.9562 1.3165 0.8553 0.6275 0.2445if (all(autovalores > 0)) {
cat("✓ S0 es positiva definida\n")
} else {
cat("✗ Reducir correlaciones a priori\n")
}## ✓ S0 es positiva definida# Grados de libertad
v0 <- p + n0 # = 15Modelo Bayesiano Multivariante
Modelo probabilístico
\[ Y_i \sim \mathcal{N}_p(\mu, \Sigma) \]
donde:
- \(\mu\) es el vector de medias
- \(\Sigma\) es la matriz de covarianza
Distribuciones a priori
Vector de medias
\[ \mu \sim \mathcal{N}_p(\mu_0, \Lambda_0) \]
Matriz de covarianza
\[ \Sigma \sim IW(v_0, S_0) \]
Objetivo científico
Los parámetros de interés son las correlaciones:
\[ \rho(T_{20}, Res_{20}) \]
y
\[ \rho(T_{150}, Res_{150}) \]
las cuales se derivan de la matriz de covarianza:
\[ \rho_{ij} = \frac{\Sigma_{ij}} {\sqrt{\Sigma_{ii}\Sigma_{jj}}} \]
cat("PARÁMETROS DE LAS DISTRIBUCIONES A PRIORI\n")## PARÁMETROS DE LAS DISTRIBUCIONES A PRIORI# -----------------------------------------------------------------------------
# 1. PRIOR DEL VECTOR DE MEDIAS mu
# -----------------------------------------------------------------------------
cat("\n")cat("PRIOR PARA EL VECTOR DE MEDIAS (mu)\n")## PRIOR PARA EL VECTOR DE MEDIAS (mu)cat("mu ~ N_p(mu0, Lambda0)\n\n")## mu ~ N_p(mu0, Lambda0)tabla_prior_mu <- data.frame(
Variable = colnames(datos_std),
Media_A_Priori = round(mu0, 3),
Varianza_A_Priori = round(diag(Lambda0), 3)
)
knitr::kable(
tabla_prior_mu,
caption = "Parámetros de la distribución Normal multivariada a priori para μ"
)| Variable | Media_A_Priori | Varianza_A_Priori |
|---|---|---|
| std_T_20cm | 0.3 | 0.025 |
| std_T_150cm | 0.2 | 0.025 |
| std_Res_20cm | -0.2 | 0.025 |
| std_Res_150cm | -0.2 | 0.025 |
| std_Altura | 0.0 | 0.025 |
# -----------------------------------------------------------------------------
# 2. MATRIZ Lambda0
# -----------------------------------------------------------------------------
cat("\n")cat("MATRIZ DE COVARIANZA PRIOR Lambda0\n")## MATRIZ DE COVARIANZA PRIOR Lambda0Lambda0_tabla <- round(as.data.frame(Lambda0), 3)
rownames(Lambda0_tabla) <- colnames(datos_std)
colnames(Lambda0_tabla) <- colnames(datos_std)
knitr::kable(
Lambda0_tabla,
caption = "Matriz de covarianza a priori Λ0"
)| std_T_20cm | std_T_150cm | std_Res_20cm | std_Res_150cm | std_Altura | |
|---|---|---|---|---|---|
| std_T_20cm | 0.025 | 0.000 | 0.000 | 0.000 | 0.000 |
| std_T_150cm | 0.000 | 0.025 | 0.000 | 0.000 | 0.000 |
| std_Res_20cm | 0.000 | 0.000 | 0.025 | 0.000 | 0.000 |
| std_Res_150cm | 0.000 | 0.000 | 0.000 | 0.025 | 0.000 |
| std_Altura | 0.000 | 0.000 | 0.000 | 0.000 | 0.025 |
# -----------------------------------------------------------------------------
# 3. PRIOR PARA LA MATRIZ DE COVARIANZA Sigma
# -----------------------------------------------------------------------------
cat("\n")cat("PRIOR PARA LA MATRIZ DE COVARIANZA (Sigma)\n")## PRIOR PARA LA MATRIZ DE COVARIANZA (Sigma)cat("Sigma ~ Inverse-Wishart(v0, S0)\n\n")## Sigma ~ Inverse-Wishart(v0, S0)tabla_prior_sigma <- data.frame(
Parametro = c(
"Grados de libertad (v0)"
),
Valor = c(v0)
)
knitr::kable(
tabla_prior_sigma,
caption = "Parámetro de grados de libertad de la prior Inverse-Wishart"
)| Parametro | Valor |
|---|---|
| Grados de libertad (v0) | 45 |
# -----------------------------------------------------------------------------
# 4. MATRIZ DE ESCALA S0
# -----------------------------------------------------------------------------
cat("\n")cat("MATRIZ DE ESCALA PRIOR S0\n")## MATRIZ DE ESCALA PRIOR S0S0_tabla <- round(as.data.frame(S0), 3)
rownames(S0_tabla) <- colnames(datos_std)
colnames(S0_tabla) <- colnames(datos_std)
knitr::kable(
S0_tabla,
caption = "Matriz de escala a priori S0 para la distribución Inverse-Wishart"
)| std_T_20cm | std_T_150cm | std_Res_20cm | std_Res_150cm | std_Altura | |
|---|---|---|---|---|---|
| std_T_20cm | 1.0 | 0.50 | -0.3 | 0.0 | -0.40 |
| std_T_150cm | 0.5 | 1.00 | 0.0 | -0.3 | -0.35 |
| std_Res_20cm | -0.3 | 0.00 | 1.0 | 0.4 | 0.00 |
| std_Res_150cm | 0.0 | -0.30 | 0.4 | 1.0 | 0.00 |
| std_Altura | -0.4 | -0.35 | 0.0 | 0.0 | 1.00 |
# -----------------------------------------------------------------------------
# 5. AUTOVALORES DE S0
# -----------------------------------------------------------------------------
tabla_autovalores <- data.frame(
Autovalor = round(eigen(S0)$values, 4)
)
knitr::kable(
tabla_autovalores,
caption = "Autovalores de S0 (verificación de positiva definida)"
)| Autovalor |
|---|
| 1.9562 |
| 1.3165 |
| 0.8553 |
| 0.6275 |
| 0.2445 |
# -----------------------------------------------------------------------------
# 7. ESTIMADORES INICIALES (ÚTILES PARA MCMC)
# -----------------------------------------------------------------------------
# Media muestral inicial
mu_init <- colMeans(datos_std)
# Covarianza muestral inicial
Sigma_init <- cov(datos_std)
cat("ESTIMADORES INICIALES\n")## ESTIMADORES INICIALEScat("Media inicial:\n")## Media inicial:print(round(mu_init, 4))## std_T_20cm std_T_150cm std_Res_20cm std_Res_150cm std_Altura
## 0 0 0 0 0cat("Covarianza inicial:\n")## Covarianza inicial:print(round(Sigma_init, 4))## std_T_20cm std_T_150cm std_Res_20cm std_Res_150cm std_Altura
## std_T_20cm 1.0000 0.8358 -0.3136 -0.1546 -0.6752
## std_T_150cm 0.8358 1.0000 -0.2712 -0.1974 -0.7399
## std_Res_20cm -0.3136 -0.2712 1.0000 0.3931 0.1821
## std_Res_150cm -0.1546 -0.1974 0.3931 1.0000 0.3605
## std_Altura -0.6752 -0.7399 0.1821 0.3605 1.0000# -----------------------------------------------------------------------------
# 8. GUARDAR OBJETOS PARA FASE 3 (MCMC / GIBBS)
# -----------------------------------------------------------------------------
save(
datos_std,
n,
p,
mu0,
Lambda0,
v0,
S0,
mu_init,
Sigma_init,
file = "objetos_modelo_bayesiano.RData"
)En esta tercera fase se realizó el cálculo analítico de las distribuciones posteriores del modelo bayesiano. Primero, se construyó la matriz de dispersión S, la cual resume la variabilidad y covariabilidad observada entre las variables mediante las sumas de cuadrados y productos cruzados respecto a la media muestral. Esta matriz constituye un componente fundamental de la verosimilitud multivariada, ya que concentra la información contenida en los 90 registros analizados. A partir de esta información y de las distribuciones a priori definidas previamente, se aplicó el teorema de Bayes para obtener la distribución posterior del vector de medias μ, modelada mediante una Normal multivariada μ∣Σ,X∼Np(μn,Λn). En esta actualización, la media posterior μn representa un equilibrio entre la información previa y la evidencia aportada por los datos, mientras que la matriz Λn refleja una reducción de la incertidumbre respecto a la prior inicial. Posteriormente, también se actualizaron los parámetros de la matriz de covarianza mediante una distribución Wishart para Σ−1, donde los grados de libertad aumentaron en función del tamaño de muestra y la matriz de escala incorporó tanto la información prior como la dispersión observada en los datos. De esta manera, las distribuciones posteriores integran formalmente el conocimiento previo con la evidencia empírica obtenida en el estudio.
# -----------------------------------------------------------------------------
# 1. CARGAR OBJETOS DE LA FASE 2
# -----------------------------------------------------------------------------
load("objetos_modelo_bayesiano.RData")
cat("FASE 3 — DISTRIBUCIONES A POSTERIORI")## FASE 3 — DISTRIBUCIONES A POSTERIORI# -----------------------------------------------------------------------------
# 2. MATRIZ DE DATOS
# -----------------------------------------------------------------------------
X <- as.matrix(datos_std)
# -----------------------------------------------------------------------------
# 3. MEDIA MUESTRAL
# -----------------------------------------------------------------------------
x_bar <- colMeans(X)
cat("MEDIA MUESTRAL\n")## MEDIA MUESTRALprint(round(x_bar, 4))## std_T_20cm std_T_150cm std_Res_20cm std_Res_150cm std_Altura
## 0 0 0 0 0# -----------------------------------------------------------------------------
# 4. MATRIZ DE DISPERSIÓN CENTRADA
# -----------------------------------------------------------------------------
# S = Σ (x_i - x_bar)(x_i - x_bar)'
# -----------------------------------------------------------------------------
S <- matrix(0, p, p)
for (i in 1:n) {
x_i <- matrix(X[i, ], ncol = 1)
x_bar_vec <- matrix(x_bar, ncol = 1)
S <- S + (x_i - x_bar_vec) %*% t(x_i - x_bar_vec)
}
cat("MATRIZ DE DISPERSIÓN S\n")## MATRIZ DE DISPERSIÓN Sprint(round(S, 4))## [,1] [,2] [,3] [,4] [,5]
## [1,] 89.0000 74.3872 -27.9073 -13.7571 -60.0920
## [2,] 74.3872 89.0000 -24.1351 -17.5686 -65.8470
## [3,] -27.9073 -24.1351 89.0000 34.9834 16.2066
## [4,] -13.7571 -17.5686 34.9834 89.0000 32.0864
## [5,] -60.0920 -65.8470 16.2066 32.0864 89.0000# -----------------------------------------------------------------------------
# 5. POSTERIOR DE mu
# -----------------------------------------------------------------------------
# p(mu | Sigma, X) ~ N_p(mu_n, Lambda_n)
#
# Lambda_n = (Lambda0^-1 + n Sigma^-1)^-1
#
# mu_n = Lambda_n (Lambda0^-1 mu0 + n Sigma^-1 x_bar)
# -----------------------------------------------------------------------------
Sigma_inv <- solve(Sigma_init)
Lambda0_inv <- solve(Lambda0)
Lambda_n <- solve(
Lambda0_inv + n * Sigma_inv
)
mu_n <- Lambda_n %*%
(
Lambda0_inv %*% matrix(mu0, ncol = 1) +
n * Sigma_inv %*% matrix(x_bar, ncol = 1)
)
# -----------------------------------------------------------------------------
# 6. POSTERIOR DE Sigma^-1
# -----------------------------------------------------------------------------
# p(Sigma^-1 | mu, X) ~ Wishart(v_n, S_n^-1)
#
# v_n = v0 + n
#
# S_n = S0 + S
# -----------------------------------------------------------------------------
v_n <- v0 + n
S_n <- S0 + S
S_n_inv <- solve(S_n)Posterior de \(\mu\)
\[ \mu \mid \Sigma, X \sim N_p(\mu_n, \Lambda_n) \]
cat("Media posterior mu_n:\n")## Media posterior mu_n:print(round(mu_n, 4))## [,1]
## std_T_20cm 0.1155
## std_T_150cm 0.1070
## std_Res_20cm -0.0972
## std_Res_150cm -0.0797
## std_Altura -0.0786cat("Covarianza posterior Lambda_n:\n")## Covarianza posterior Lambda_n:print(round(Lambda_n, 4))## std_T_20cm std_T_150cm std_Res_20cm std_Res_150cm std_Altura
## std_T_20cm 0.0059 0.0041 -0.0013 -0.0002 -0.0029
## std_T_150cm 0.0041 0.0057 -0.0010 -0.0005 -0.0034
## std_Res_20cm -0.0013 -0.0010 0.0072 0.0020 0.0003
## std_Res_150cm -0.0002 -0.0005 0.0020 0.0072 0.0017
## std_Altura -0.0029 -0.0034 0.0003 0.0017 0.0061Posterior de \(\Sigma^{-1}\)
\[ \Sigma^{-1} \mid \mu, X \sim \text{Wishart}(v_n, S_n^{-1}) \]
cat("Grados de libertad posteriores v_n:\n")## Grados de libertad posteriores v_n:print(v_n)## [1] 135cat("Matriz S_n:\n")## Matriz S_n:print(round(S_n, 4))## std_T_20cm std_T_150cm std_Res_20cm std_Res_150cm std_Altura
## std_T_20cm 90.0000 74.8872 -28.2073 -13.7571 -60.4920
## std_T_150cm 74.8872 90.0000 -24.1351 -17.8686 -66.1970
## std_Res_20cm -28.2073 -24.1351 90.0000 35.3834 16.2066
## std_Res_150cm -13.7571 -17.8686 35.3834 90.0000 32.0864
## std_Altura -60.4920 -66.1970 16.2066 32.0864 90.0000# -----------------------------------------------------------------------------
# 7. MATRIZ DE CORRELACIÓN POSTERIOR INICIAL
# -----------------------------------------------------------------------------
# Conversión:
#
# rho_ij = Sigma_ij / sqrt(Sigma_ii * Sigma_jj)
# -----------------------------------------------------------------------------
Sigma_post_init <- S_n / v_n
R_post_init <- cov2cor(Sigma_post_init)
cat("MATRIZ DE CORRELACIÓN POSTERIOR INICIAL\n")## MATRIZ DE CORRELACIÓN POSTERIOR INICIALprint(round(R_post_init, 4))## std_T_20cm std_T_150cm std_Res_20cm std_Res_150cm std_Altura
## std_T_20cm 1.0000 0.8321 -0.3134 -0.1529 -0.6721
## std_T_150cm 0.8321 1.0000 -0.2682 -0.1985 -0.7355
## std_Res_20cm -0.3134 -0.2682 1.0000 0.3931 0.1801
## std_Res_150cm -0.1529 -0.1985 0.3931 1.0000 0.3565
## std_Altura -0.6721 -0.7355 0.1801 0.3565 1.0000# -----------------------------------------------------------------------------
# 8. CORRELACIONES CIENTÍFICAS DE INTERÉS
# -----------------------------------------------------------------------------
rho_post_T20_Res20 <-
R_post_init[idx_T20, idx_Res20]
rho_post_T150_Res150 <-
R_post_init[idx_T150, idx_Res150]
# Tabla de correlaciones posteriores
tabla_corr <- data.frame(
Variable = c("rho(T20, Res20)",
"rho(T150, Res150)"),
Valor = c(round(rho_post_T20_Res20, 4),
round(rho_post_T150_Res150, 4))
)
knitr::kable(tabla_corr,
caption = "Correlaciones posteriores iniciales")| Variable | Valor |
|---|---|
| rho(T20, Res20) | -0.3134 |
| rho(T150, Res150) | -0.1985 |
Resumen del Modelo Posterior
Posterior de \(\mu\)
\[ p(\mu \mid \Sigma, X) \sim N_p(\mu_n, \Lambda_n) \]
Posterior de \(\Sigma^{-1}\)
\[ p(\Sigma^{-1} \mid \mu, X) \sim \text{Wishart}(v_n, S_n^{-1}) \]
Objetivo científico
Inferir las distribuciones posteriores de las correlaciones:
\[ \rho (T20, Res20) \]
\[ \rho(T150, Res150) \]
# -----------------------------------------------------------------------------
# 10. GUARDAR OBJETOS PARA FASE 4 (GIBBS SAMPLER)
# -----------------------------------------------------------------------------
save(
# Datos
X,
x_bar,
S,
# Posteriores
mu_n,
Lambda_n,
v_n,
S_n,
S_n_inv,
# Correlaciones
R_post_init,
rho_post_T20_Res20,
rho_post_T150_Res150,
# Índices
idx_T20,
idx_Res20,
idx_T150,
idx_Res150,
file = "objetos_posterior_bayesiana.RData"
)En esta cuarta fase se implementó el algoritmo de Gibbs Sampling mediante técnicas MCMC para generar muestras de las distribuciones posteriores y realizar la inferencia bayesiana completa. Aunque las distribuciones posteriores de los parámetros pueden obtenerse analíticamente, el muestreo iterativo permite explorar su comportamiento y estimar cantidades de interés a partir de simulaciones. Para ello, se ejecutaron 10.000 iteraciones del algoritmo, descartando las primeras 2.000 como período de burn-in, ya que durante esta etapa la cadena aún no alcanza estabilidad en la región de mayor probabilidad posterior. En cada iteración, el algoritmo alternó dos pasos condicionales: primero, se muestreó el vector de medias μ condicionado a la matriz de covarianza actual Σ y a los datos, utilizando una distribución Normal multivariada; posteriormente, se muestreó Σ condicionado al nuevo valor de μ, actualizando la matriz de dispersión y generando muestras a partir de una distribución Wishart. A partir de cada matriz de covarianza simulada se calculó la correspondiente matriz de correlación, extrayendo específicamente las correlaciones entre temperatura superficial y resistividad superficial, así como entre temperatura profunda y resistividad profunda. Finalmente, tras eliminar las iteraciones iniciales de estabilización, se conservaron 8.000 muestras posteriores representativas, las cuales constituyen la base para toda la inferencia estadística y el análisis de incertidumbre del modelo.
# -----------------------------------------------------------------------------
# 1. CARGAR OBJETOS DE LA FASE 3
# -----------------------------------------------------------------------------
load("objetos_posterior_bayesiana.RData")
load("objetos_modelo_bayesiano.RData")
cat("FASE 4 — GIBBS SAMPLER (MCMC)\n")## FASE 4 — GIBBS SAMPLER (MCMC)# -----------------------------------------------------------------------------
# 2. CONFIGURACIÓN MCMC
# -----------------------------------------------------------------------------
n_iter <- 100000
burn_in <- 2000
cat("Número total de iteraciones:", n_iter, "\n")## Número total de iteraciones: 1e+05cat("Burn-in:", burn_in, "\n")## Burn-in: 2000# -----------------------------------------------------------------------------
# 3. OBJETOS PARA ALMACENAR MUESTRAS
# -----------------------------------------------------------------------------
# Muestras del vector de medias
mu_samples <- array(NA, dim = c(n_iter, p))
# Muestras de matrices de covarianza
Sigma_samples <- array(NA, dim = c(p, p, n_iter))
# Correlaciones de interés
rho_T20_Res20 <- numeric(n_iter)
rho_T150_Res150 <- numeric(n_iter)
# -----------------------------------------------------------------------------
# 4. INICIALIZACIÓN
# -----------------------------------------------------------------------------
mu_current <- mu_init
Sigma_current <- Sigma_init
# -----------------------------------------------------------------------------
# 5. GIBBS SAMPLER
# -----------------------------------------------------------------------------
for (iter in 1:n_iter) {
# ===========================================================================
# PASO 1: MUESTREAR mu | Sigma, X
# ===========================================================================
Sigma_inv <- solve(Sigma_current)
Lambda_n <- solve(
solve(Lambda0) + n * Sigma_inv
)
mu_n <- Lambda_n %*%
(
solve(Lambda0) %*% matrix(mu0, ncol = 1) +
n * Sigma_inv %*% matrix(x_bar, ncol = 1)
)
mu_current <- as.vector(
rmvnorm(
n = 1,
mean = as.vector(mu_n),
sigma = Lambda_n
)
)
# ===========================================================================
# PASO 2: MUESTREAR Sigma | mu, X
# ===========================================================================
S_mu <- matrix(0, p, p)
for (i in 1:n) {
x_i <- matrix(X[i, ], ncol = 1)
mu_vec <- matrix(mu_current, ncol = 1)
S_mu <- S_mu + (x_i - mu_vec) %*% t(x_i - mu_vec)
}
S_n <- S0 + S_mu
v_n <- v0 + n
# Muestrear Wishart
W_sample <- rwish(
v = v_n,
S = solve(S_n)
)
# Convertir a Inverse-Wishart
Sigma_current <- solve(W_sample)
# ===========================================================================
# GUARDAR MUESTRAS
# ===========================================================================
mu_samples[iter, ] <- mu_current
Sigma_samples[, , iter] <- Sigma_current
# ===========================================================================
# CALCULAR CORRELACIONES POSTERIORES
# ===========================================================================
R_current <- cov2cor(Sigma_current)
rho_T20_Res20[iter] <-
R_current[idx_T20, idx_Res20]
rho_T150_Res150[iter] <-
R_current[idx_T150, idx_Res150]
# ===========================================================================
# PROGRESO
# ===========================================================================
if (iter %% 1000 == 0) {
cat("Iteración:", iter, "de", n_iter, "\n")
}
}## Iteración: 1000 de 1e+05
## Iteración: 2000 de 1e+05
## Iteración: 3000 de 1e+05
## Iteración: 4000 de 1e+05
## Iteración: 5000 de 1e+05
## Iteración: 6000 de 1e+05
## Iteración: 7000 de 1e+05
## Iteración: 8000 de 1e+05
## Iteración: 9000 de 1e+05
## Iteración: 10000 de 1e+05
## Iteración: 11000 de 1e+05
## Iteración: 12000 de 1e+05
## Iteración: 13000 de 1e+05
## Iteración: 14000 de 1e+05
## Iteración: 15000 de 1e+05
## Iteración: 16000 de 1e+05
## Iteración: 17000 de 1e+05
## Iteración: 18000 de 1e+05
## Iteración: 19000 de 1e+05
## Iteración: 20000 de 1e+05
## Iteración: 21000 de 1e+05
## Iteración: 22000 de 1e+05
## Iteración: 23000 de 1e+05
## Iteración: 24000 de 1e+05
## Iteración: 25000 de 1e+05
## Iteración: 26000 de 1e+05
## Iteración: 27000 de 1e+05
## Iteración: 28000 de 1e+05
## Iteración: 29000 de 1e+05
## Iteración: 30000 de 1e+05
## Iteración: 31000 de 1e+05
## Iteración: 32000 de 1e+05
## Iteración: 33000 de 1e+05
## Iteración: 34000 de 1e+05
## Iteración: 35000 de 1e+05
## Iteración: 36000 de 1e+05
## Iteración: 37000 de 1e+05
## Iteración: 38000 de 1e+05
## Iteración: 39000 de 1e+05
## Iteración: 40000 de 1e+05
## Iteración: 41000 de 1e+05
## Iteración: 42000 de 1e+05
## Iteración: 43000 de 1e+05
## Iteración: 44000 de 1e+05
## Iteración: 45000 de 1e+05
## Iteración: 46000 de 1e+05
## Iteración: 47000 de 1e+05
## Iteración: 48000 de 1e+05
## Iteración: 49000 de 1e+05
## Iteración: 50000 de 1e+05
## Iteración: 51000 de 1e+05
## Iteración: 52000 de 1e+05
## Iteración: 53000 de 1e+05
## Iteración: 54000 de 1e+05
## Iteración: 55000 de 1e+05
## Iteración: 56000 de 1e+05
## Iteración: 57000 de 1e+05
## Iteración: 58000 de 1e+05
## Iteración: 59000 de 1e+05
## Iteración: 60000 de 1e+05
## Iteración: 61000 de 1e+05
## Iteración: 62000 de 1e+05
## Iteración: 63000 de 1e+05
## Iteración: 64000 de 1e+05
## Iteración: 65000 de 1e+05
## Iteración: 66000 de 1e+05
## Iteración: 67000 de 1e+05
## Iteración: 68000 de 1e+05
## Iteración: 69000 de 1e+05
## Iteración: 70000 de 1e+05
## Iteración: 71000 de 1e+05
## Iteración: 72000 de 1e+05
## Iteración: 73000 de 1e+05
## Iteración: 74000 de 1e+05
## Iteración: 75000 de 1e+05
## Iteración: 76000 de 1e+05
## Iteración: 77000 de 1e+05
## Iteración: 78000 de 1e+05
## Iteración: 79000 de 1e+05
## Iteración: 80000 de 1e+05
## Iteración: 81000 de 1e+05
## Iteración: 82000 de 1e+05
## Iteración: 83000 de 1e+05
## Iteración: 84000 de 1e+05
## Iteración: 85000 de 1e+05
## Iteración: 86000 de 1e+05
## Iteración: 87000 de 1e+05
## Iteración: 88000 de 1e+05
## Iteración: 89000 de 1e+05
## Iteración: 90000 de 1e+05
## Iteración: 91000 de 1e+05
## Iteración: 92000 de 1e+05
## Iteración: 93000 de 1e+05
## Iteración: 94000 de 1e+05
## Iteración: 95000 de 1e+05
## Iteración: 96000 de 1e+05
## Iteración: 97000 de 1e+05
## Iteración: 98000 de 1e+05
## Iteración: 99000 de 1e+05
## Iteración: 100000 de 1e+05# -----------------------------------------------------------------------------
# 6. ELIMINAR BURN-IN
# -----------------------------------------------------------------------------
mu_post <- mu_samples[(burn_in + 1):n_iter, ]
Sigma_post <- Sigma_samples[, , (burn_in + 1):n_iter]
rho20_post <- rho_T20_Res20[(burn_in + 1):n_iter]
rho150_post <- rho_T150_Res150[(burn_in + 1):n_iter]
cat("Número final de muestras posteriores:",
length(rho20_post), "\n")## Número final de muestras posteriores: 98000# -----------------------------------------------------------------------------
# 7. RESUMEN POSTERIOR DE CORRELACIONES
# -----------------------------------------------------------------------------
library(knitr)
cat("\n")cat("RESUMEN POSTERIOR — CORRELACIONES\n")## RESUMEN POSTERIOR — CORRELACIONES# -----------------------------------------------------------------------------
# TABLA 1 — T20 vs Res20
# -----------------------------------------------------------------------------
tabla_rho20 <- data.frame(
Estadistico = c(
"Media posterior",
"Desviación estándar",
"Percentil 2.5%",
"Mediana",
"Percentil 97.5%",
"% correlaciones negativas"
),
Valor = c(
mean(rho20_post),
sd(rho20_post),
quantile(rho20_post, 0.025),
median(rho20_post),
quantile(rho20_post, 0.975),
mean(rho20_post < 0) * 100
)
)
tabla_rho20$Valor <- round(tabla_rho20$Valor, 4)
cat("\n")cat("Correlación posterior: T20 vs Res20\n")## Correlación posterior: T20 vs Res20kable(
tabla_rho20,
caption = "Resumen posterior de la correlación entre temperatura y resistividad a 20 cm"
)| Estadistico | Valor |
|---|---|
| Media posterior | -0.3167 |
| Desviación estándar | 0.0786 |
| Percentil 2.5% | -0.4649 |
| Mediana | -0.3191 |
| Percentil 97.5% | -0.1568 |
| % correlaciones negativas | 99.9878 |
# -----------------------------------------------------------------------------
# TABLA 2 — T150 vs Res150
# -----------------------------------------------------------------------------
tabla_rho150 <- data.frame(
Estadistico = c(
"Media posterior",
"Desviación estándar",
"Percentil 2.5%",
"Mediana",
"Percentil 97.5%",
"% correlaciones negativas"
),
Valor = c(
mean(rho150_post),
sd(rho150_post),
quantile(rho150_post, 0.025),
median(rho150_post),
quantile(rho150_post, 0.975),
mean(rho150_post < 0) * 100
)
)
tabla_rho150$Valor <- round(tabla_rho150$Valor, 4)
cat("\n")cat("Correlación posterior: T150 vs Res150\n")## Correlación posterior: T150 vs Res150kable(
tabla_rho150,
caption = "Resumen posterior de la correlación entre temperatura y resistividad a 150 cm"
)| Estadistico | Valor |
|---|---|
| Media posterior | -0.2017 |
| Desviación estándar | 0.0838 |
| Percentil 2.5% | -0.3614 |
| Mediana | -0.2029 |
| Percentil 97.5% | -0.0324 |
| % correlaciones negativas | 99.0408 |
# -----------------------------------------------------------------------------
# 8. GUARDAR RESULTADOS
# -----------------------------------------------------------------------------
save(
mu_post,
Sigma_post,
rho20_post,
rho150_post,
tabla_rho20,
tabla_rho150,
file = "resultados_mcmc_gibbs.RData"
)En esta quinta fase se evaluó la convergencia y calidad de las cadenas generadas por el algoritmo MCMC, paso fundamental antes de interpretar los resultados posteriores. Para ello, se utilizaron diferentes diagnósticos de convergencia. En primer lugar, los traceplots permitieron visualizar la evolución de las cadenas a lo largo de las iteraciones; una cadena adecuadamente convergida debe comportarse como una fluctuación aleatoria alrededor de un valor estable, sin presentar tendencias sistemáticas. Posteriormente, se analizó la autocorrelación entre muestras consecutivas, con el fin de determinar el grado de dependencia dentro de la cadena. Valores bajos de autocorrelación indican un muestreo más eficiente y cercano a la independencia entre iteraciones, mientras que valores altos sugieren una exploración lenta de la distribución posterior. Finalmente, se calculó el tamaño de muestra efectivo (Effective Sample Size, ESS), que estima cuántas observaciones independientes equivalen realmente las muestras obtenidas después del burn-in. Un ESS elevado indica que la cadena proporciona estimaciones estables y confiables, mientras que valores bajos sugieren la necesidad de aumentar el número de iteraciones para mejorar la precisión de la inferencia bayesiana.
# -----------------------------------------------------------------------------
# 1. CARGAR RESULTADOS MCMC
# -----------------------------------------------------------------------------
load("resultados_mcmc_gibbs.RData")
cat("FASE 5 — DIAGNÓSTICO DE CONVERGENCIA\n")## FASE 5 — DIAGNÓSTICO DE CONVERGENCIA# -----------------------------------------------------------------------------
# 2. CONVERTIR A OBJETOS MCMC
# -----------------------------------------------------------------------------
mcmc_rho20 <- mcmc(rho20_post)
mcmc_rho150 <- mcmc(rho150_post)
# -----------------------------------------------------------------------------
# 3. TRACEPLOTS
# -----------------------------------------------------------------------------
df_rho20 <- data.frame(
Iteracion = 1:length(rho20_post),
rho = rho20_post
)
df_rho150 <- data.frame(
Iteracion = 1:length(rho150_post),
rho = rho150_post
)
# -------------------------------------------------------------------------
# Traceplot rho(T20, Res20)
# -------------------------------------------------------------------------
p1 <- ggplot(df_rho20,
aes(x = Iteracion, y = rho)) +
geom_line(color = "#2E86AB", linewidth = 0.3) +
labs(
title = expression(paste("Traceplot: ", rho,
"(T20, Res20)")),
x = "Iteración",
y = expression(rho)
) +
theme_minimal()
# -------------------------------------------------------------------------
# Traceplot rho(T150, Res150)
# -------------------------------------------------------------------------
p2 <- ggplot(df_rho150,
aes(x = Iteracion, y = rho)) +
geom_line(color = "#E76F51", linewidth = 0.3) +
labs(
title = expression(paste("Traceplot: ", rho,
"(T150, Res150)")),
x = "Iteración",
y = expression(rho)
) +
theme_minimal()
# -------------------------------------------------------------------------
# Guardar figura
# -------------------------------------------------------------------------
fig_trace <- grid.arrange(p1, p2, ncol = 1)
ggsave(
"fig_traceplots_mcmc.png",
fig_trace,
width = 10,
height = 7,
dpi = 180
)
# -----------------------------------------------------------------------------
# 4. AUTOCORRELACIÓN
# -----------------------------------------------------------------------------
library(knitr)
acf_rho20 <- autocorr.diag(mcmc_rho20)
acf_rho150 <- autocorr.diag(mcmc_rho150)
# Crear tabla de autocorrelaciones
tabla_acf <- data.frame(
Lag = 1:10,
`rho(T20, Res20)` = round(acf_rho20[1:5], 4),
`rho(T150, Res150)` = round(acf_rho150[1:5], 4)
)
cat("\n")cat("AUTOCORRELACIÓN MCMC\n")## AUTOCORRELACIÓN MCMCkable(
tabla_acf,
caption = "Autocorrelaciones de las cadenas MCMC para las correlaciones posteriores"
)| Lag | rho.T20..Res20. | rho.T150..Res150. |
|---|---|---|
| 1 | 1.0000 | 1.0000 |
| 2 | 0.0079 | 0.0076 |
| 3 | -0.0002 | 0.0006 |
| 4 | -0.0061 | 0.0010 |
| 5 | -0.0036 | 0.0028 |
| 6 | 1.0000 | 1.0000 |
| 7 | 0.0079 | 0.0076 |
| 8 | -0.0002 | 0.0006 |
| 9 | -0.0061 | 0.0010 |
| 10 | -0.0036 | 0.0028 |
# -----------------------------------------------------------------------------
# 5. EFFECTIVE SAMPLE SIZE (ESS)
# -----------------------------------------------------------------------------
ess_rho20 <- effectiveSize(mcmc_rho20)
ess_rho150 <- effectiveSize(mcmc_rho150)
tabla_ess <- data.frame(
Parametro = c(
"rho(T20, Res20)",
"rho(T150, Res150)"
),
ESS = c(
round(ess_rho20, 2),
round(ess_rho150, 2)
)
)
cat("\n")cat("EFFECTIVE SAMPLE SIZE (ESS)\n")## EFFECTIVE SAMPLE SIZE (ESS)kable(
tabla_ess,
caption = "Tamaño efectivo de muestra de las cadenas MCMC"
)| Parametro | ESS |
|---|---|
| rho(T20, Res20) | 97487.25 |
| rho(T150, Res150) | 97596.94 |
# -----------------------------------------------------------------------------
# 6. INTERPRETACIÓN AUTOMÁTICA
# -----------------------------------------------------------------------------
interpretacion_rho20 <- ifelse(
ess_rho20 > 400,
"ESS adecuado (> 400)",
"ESS bajo — aumentar iteraciones"
)
interpretacion_rho150 <- ifelse(
ess_rho150 > 400,
"ESS adecuado (> 400)",
"ESS bajo — aumentar iteraciones"
)
tabla_diag <- data.frame(
Parametro = c(
"rho(T20, Res20)",
"rho(T150, Res150)"
),
ESS = c(
round(ess_rho20, 2),
round(ess_rho150, 2)
),
Diagnostico = c(
interpretacion_rho20,
interpretacion_rho150
)
)
cat("\n")cat("INTERPRETACIÓN DEL DIAGNÓSTICO\n")## INTERPRETACIÓN DEL DIAGNÓSTICOkable(
tabla_diag,
caption = "Diagnóstico automático de convergencia MCMC"
)| Parametro | ESS | Diagnostico |
|---|---|---|
| rho(T20, Res20) | 97487.25 | ESS adecuado (> 400) |
| rho(T150, Res150) | 97596.94 | ESS adecuado (> 400) |
# -----------------------------------------------------------------------------
# 7. GUARDAR OBJETOS DE DIAGNÓSTICO
# -----------------------------------------------------------------------------
save(
mcmc_rho20,
mcmc_rho150,
ess_rho20,
ess_rho150,
acf_rho20,
acf_rho150,
file = "diagnostico_convergencia.RData"
)En esta sexta fase se realizó la inferencia bayesiana final y la interpretación de los resultados obtenidos a partir de las muestras posteriores generadas por el Gibbs Sampler. Primero, se calcularon intervalos de credibilidad al 95% para cada componente del vector de medias μ utilizando los percentiles 2.5% y 97.5% de las muestras posteriores, proporcionando rangos probabilísticos donde es más probable que se encuentren los valores reales de los parámetros. Posteriormente, se estimó la probabilidad posterior de correlaciones negativas entre variables de interés, calculando la proporción de muestras posteriores menores que cero; este procedimiento permite realizar afirmaciones probabilísticas directas sobre la dirección de la relación entre temperatura y resistividad, una ventaja distintiva del enfoque bayesiano frente al clásico. Además, se compararon las correlaciones posteriores con las estimaciones clásicas obtenidas mediante máxima verosimilitud (MLE), observándose generalmente resultados similares debido al tamaño de muestra y al uso de priors poco informativas, aunque el enfoque bayesiano aporta adicionalmente una cuantificación explícita de la incertidumbre. Finalmente, se construyó la distribución predictiva posterior generando observaciones sintéticas de nuevos puntos de campo a partir de muestras de (μ,Σ) obtenidas de la posterior. Esto permitió incorporar la incertidumbre de los parámetros y estimar qué valores futuros de temperatura y resistividad podrían esperarse en la zona de estudio junto con su variabilidad asociada.
# -----------------------------------------------------------------------------
# 1. CARGAR RESULTADOS
# -----------------------------------------------------------------------------
load("resultados_mcmc_gibbs.RData")
load("objetos_modelo_bayesiano.RData")
# -----------------------------------------------------------------------------
# FASE 6 — INFERENCIA Y RESULTADOS
# -----------------------------------------------------------------------------
variables <- colnames(datos_std)
cat("FASE 6 — INFERENCIA Y RESULTADOS\n")## FASE 6 — INFERENCIA Y RESULTADOS# -----------------------------------------------------------------------------
# 1. RECONSTRUIR MEDIAS Y SDs ORIGINALES DE LAS VARIABLES TRANSFORMADAS
# (los mismos valores que usó scale() internamente)
# -----------------------------------------------------------------------------
# Recrear datos_transf igual que en Fase 1
datos_transf_reconstruido <- datos # datos originales cargados en Fase 1
for (v in variables_originales) {
if (resultados_diag[[v]]$necesita_log) {
nuevo_nombre <- paste0("log10_", v)
datos_transf_reconstruido[[nuevo_nombre]] <- log10(datos_transf_reconstruido[[v]])
}
}
# cols_analisis: las columnas que entraron a scale()
cols_analisis_rec <- sapply(variables_originales, function(v) {
if (resultados_diag[[v]]$necesita_log) paste0("log10_", v) else v
})
# Medias y SDs que usó scale() — en escala transformada (log10 donde aplique)
medias_transf <- colMeans(datos_transf_reconstruido[, cols_analisis_rec])
sds_transf <- apply(datos_transf_reconstruido[, cols_analisis_rec], 2, sd)
cat("Medias usadas en la estandarización (escala transformada):\n")## Medias usadas en la estandarización (escala transformada):print(round(medias_transf, 6))## log10_T_20cm T_150cm log10_Res_20cm log10_Res_150cm log10_Altura
## 1.171281 14.862222 2.281559 2.233800 3.485080cat("SDs usadas en la estandarización:\n")## SDs usadas en la estandarización:print(round(sds_transf, 6))## log10_T_20cm T_150cm log10_Res_20cm log10_Res_150cm log10_Altura
## 0.045852 1.376331 0.410142 0.526703 0.010950# -----------------------------------------------------------------------------
# 2. FUNCIÓN DE DESESTANDARIZACIÓN
# Deshace: z = (x_transf - media_transf) / sd_transf
# Luego: si la variable tenía log10 → x_original = 10^x_transf
# -----------------------------------------------------------------------------
desestandarizar_mu_post <- function(mu_post_std,
medias_transf,
sds_transf,
variables_originales,
resultados_diag) {
n_muestras <- nrow(mu_post_std)
n_vars <- ncol(mu_post_std)
mu_post_orig <- matrix(NA, nrow = n_muestras, ncol = n_vars)
colnames(mu_post_orig) <- variables_originales
for (j in seq_len(n_vars)) {
v <- variables_originales[j]
# Paso 1 — deshacer la estandarización (volver a escala log10 o escala original)
mu_transf <- mu_post_std[, j] * sds_transf[j] + medias_transf[j]
# Paso 2 — deshacer log10 si aplica
if (resultados_diag[[v]]$necesita_log) {
mu_post_orig[, j] <- 10^mu_transf
} else {
mu_post_orig[, j] <- mu_transf
}
}
return(mu_post_orig)
}
# -----------------------------------------------------------------------------
# 3. APLICAR DESESTANDARIZACIÓN A mu_post
# -----------------------------------------------------------------------------
mu_post_original <- desestandarizar_mu_post(
mu_post_std = mu_post,
medias_transf = medias_transf,
sds_transf = sds_transf,
variables_originales = variables_originales,
resultados_diag = resultados_diag
)
# -----------------------------------------------------------------------------
# 4. RESUMEN POSTERIOR EN ESCALA ORIGINAL
# -----------------------------------------------------------------------------
media_post_mu_orig <- apply(mu_post_original, 2, mean)
ic95_mu_orig <- apply(
mu_post_original,
2,
quantile,
probs = c(0.025, 0.975)
)
# Unidades reales de cada variable
unidades <- c(
T_20cm = "°C",
T_150cm = "°C",
Res_20cm = "Ω·m",
Res_150cm = "Ω·m",
Altura = "m"
)
tabla_mu_orig <- data.frame(
Variable = variables_originales,
Unidad = unidades[variables_originales],
Media_Posterior = round(media_post_mu_orig, 3),
IC95_Lower = round(ic95_mu_orig[1, ], 3),
IC95_Upper = round(ic95_mu_orig[2, ], 3)
)
cat("\n")cat("Distribuciones posteriores de mu — ESCALA ORIGINAL\n")## Distribuciones posteriores de mu — ESCALA ORIGINALkable(
tabla_mu_orig,
caption = "Media posterior e intervalos de credibilidad al 95% para el vector de medias (escala original)"
)| Variable | Unidad | Media_Posterior | IC95_Lower | IC95_Upper | |
|---|---|---|---|---|---|
| T_20cm | T_20cm | °C | 14.986 | 14.772 | 15.204 |
| T_150cm | T_150cm | °C | 14.986 | 14.802 | 15.173 |
| Res_20cm | Res_20cm | Ω·m | 177.827 | 154.191 | 203.766 |
| Res_150cm | Res_150cm | Ω·m | 158.861 | 132.080 | 189.174 |
| Altura | Altura | m | 3050.275 | 3039.695 | 3060.821 |
# -----------------------------------------------------------------------------
# 2. CORRELACIONES POSTERIORES
# -----------------------------------------------------------------------------
media_rho20 <- mean(rho20_post)
ic95_rho20 <- quantile(
rho20_post,
c(0.025, 0.975)
)
media_rho150 <- mean(rho150_post)
ic95_rho150 <- quantile(
rho150_post,
c(0.025, 0.975)
)
tabla_corr <- data.frame(
Parametro = c(
"rho(T20, Res20)",
"rho(T150, Res150)"
),
Media_Posterior = round(
c(media_rho20, media_rho150),
4
),
IC95_Lower = round(
c(ic95_rho20[1], ic95_rho150[1]),
4
),
IC95_Upper = round(
c(ic95_rho20[2], ic95_rho150[2]),
4
)
)
cat("\n")cat("Correlaciones posteriores\n")## Correlaciones posterioreskable(
tabla_corr,
caption = "Medias posteriores e intervalos de credibilidad para las correlaciones geofísicas"
)| Parametro | Media_Posterior | IC95_Lower | IC95_Upper |
|---|---|---|---|
| rho(T20, Res20) | -0.3167 | -0.4649 | -0.1568 |
| rho(T150, Res150) | -0.2017 | -0.3614 | -0.0324 |
# -----------------------------------------------------------------------------
# 3. PROBABILIDAD POSTERIOR DE CORRELACIÓN NEGATIVA
# -----------------------------------------------------------------------------
prob_neg_rho20 <- mean(rho20_post < 0)
prob_neg_rho150 <- mean(rho150_post < 0)
tabla_prob <- data.frame(
Parametro = c(
"rho(T20, Res20)",
"rho(T150, Res150)"
),
Probabilidad_Negativa = round(
c(prob_neg_rho20, prob_neg_rho150),
4
)
)
cat("\n")cat("Probabilidad posterior de correlación negativa\n")## Probabilidad posterior de correlación negativakable(
tabla_prob,
caption = "Probabilidad posterior de correlación negativa"
)| Parametro | Probabilidad_Negativa |
|---|---|
| rho(T20, Res20) | 0.9999 |
| rho(T150, Res150) | 0.9904 |
# -----------------------------------------------------------------------------
# 4. COMPARACIÓN BAYES VS MLE CLÁSICO
# -----------------------------------------------------------------------------
R_mle <- cor(datos_std)
rho20_mle <- R_mle[idx_T20, idx_Res20]
rho150_mle <- R_mle[idx_T150, idx_Res150]
tabla_comp <- data.frame(
Parametro = c(
"rho(T20, Res20)",
"rho(T150, Res150)"
),
MLE_Clasico = c(
rho20_mle,
rho150_mle
),
Media_Posterior_Bayes = c(
media_rho20,
media_rho150
)
)
tabla_comp[, -1] <- round(tabla_comp[, -1], 4)
cat("\n")cat("Comparación: Bayesiano vs MLE clásico\n")## Comparación: Bayesiano vs MLE clásicokable(
tabla_comp,
caption = "Comparación entre estimaciones clásicas y bayesianas"
)| Parametro | MLE_Clasico | Media_Posterior_Bayes |
|---|---|---|
| rho(T20, Res20) | -0.3136 | -0.3167 |
| rho(T150, Res150) | -0.1974 | -0.2017 |
# -----------------------------------------------------------------------------
# 5. DISTRIBUCIÓN PREDICTIVA POSTERIOR
# -----------------------------------------------------------------------------
n_pred <- 1000
Y_pred <- matrix(NA, nrow = n_pred, ncol = p)
for (i in 1:n_pred) {
idx <- sample(1:dim(Sigma_post)[3], 1)
mu_i <- mu_post[idx, ]
Sigma_i <- Sigma_post[, , idx]
Y_pred[i, ] <- rmvnorm(
n = 1,
mean = mu_i,
sigma = Sigma_i
)
}
colnames(Y_pred) <- variables
# -----------------------------------------------------------------------------
# 6. RESUMEN PREDICTIVO
# -----------------------------------------------------------------------------
tabla_pred <- data.frame(
Variable = variables,
Media_Predictiva =
round(colMeans(Y_pred), 4),
SD_Predictiva =
round(apply(Y_pred, 2, sd), 4)
)
cat("\n")cat("Distribución predictiva posterior\n")## Distribución predictiva posteriorkable(
tabla_pred,
caption = "Resumen de la distribución predictiva posterior"
)| Variable | Media_Predictiva | SD_Predictiva | |
|---|---|---|---|
| std_T_20cm | std_T_20cm | 0.0732 | 0.8653 |
| std_T_150cm | std_T_150cm | 0.0662 | 0.8449 |
| std_Res_20cm | std_Res_20cm | -0.0724 | 0.8685 |
| std_Res_150cm | std_Res_150cm | -0.0857 | 0.8635 |
| std_Altura | std_Altura | -0.0658 | 0.8333 |
# -----------------------------------------------------------------------------
# 8. FIGURAS POSTERIORES
# -----------------------------------------------------------------------------
# -------------------------------------------------------------------------
# Histograma rho20
# -------------------------------------------------------------------------
df_rho20 <- data.frame(rho = rho20_post)
p1 <- ggplot(df_rho20,
aes(x = rho)) +
geom_histogram(
bins = 35,
fill = "#2E86AB",
color = "white",
alpha = 0.8
) +
geom_vline(
xintercept = mean(rho20_post),
color = "red",
linewidth = 1
) +
labs(
title = expression(paste(
"Posterior: ", rho,
"(T20, Res20)"
)),
x = expression(rho),
y = "Frecuencia"
) +
theme_minimal()
# -------------------------------------------------------------------------
# Histograma rho150
# -------------------------------------------------------------------------
df_rho150 <- data.frame(rho = rho150_post)
p2 <- ggplot(df_rho150,
aes(x = rho)) +
geom_histogram(
bins = 35,
fill = "#E76F51",
color = "white",
alpha = 0.8
) +
geom_vline(
xintercept = mean(rho150_post),
color = "red",
linewidth = 1
) +
labs(
title = expression(paste(
"Posterior: ", rho,
"(T150, Res150)"
)),
x = expression(rho),
y = "Frecuencia"
) +
theme_minimal()
# -------------------------------------------------------------------------
# Guardar figura
# -------------------------------------------------------------------------
fig_post <- grid.arrange(
p1,
p2,
ncol = 1
)
ggsave(
"fig_posteriores_correlacion.png",
fig_post,
width = 9,
height = 7,
dpi = 180
)
# -----------------------------------------------------------------------------
# 9. GUARDAR RESULTADOS FINALES
# -----------------------------------------------------------------------------
save(
tabla_mu_orig, # <-- nombre actualizado
tabla_comp,
tabla_pred,
media_rho20,
media_rho150,
ic95_rho20,
ic95_rho150,
prob_neg_rho20,
prob_neg_rho150,
Y_pred,
file = "resultados_finales_bayesianos.RData"
)
# -----------------------------------------------------------------------------
# FIGURA — COMPARACIÓN PRIOR vs POSTERIOR
# para rho(T20, Res20) y rho(T150, Res150)
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# 1. SIMULAR MUESTRAS DE LA DISTRIBUCIÓN PRIOR DE LAS CORRELACIONES
# -----------------------------------------------------------------------------
# La prior de Sigma es Inverse-Wishart(v0, S0)
# De ahí se extraen correlaciones prior simulando igual que en el Gibbs Sampler
# pero SIN actualizar con los datos — solo usando v0 y S0
# -----------------------------------------------------------------------------
n_prior_samples <- 98000
rho20_prior <- numeric(n_prior_samples)
rho150_prior <- numeric(n_prior_samples)
for (i in 1:n_prior_samples) {
# Muestrear Sigma de la prior Inverse-Wishart(v0, S0)
W_prior <- rwish(
v = v0,
S = solve(S0)
)
Sigma_prior_i <- solve(W_prior)
# Convertir a correlación
R_prior_i <- cov2cor(Sigma_prior_i)
rho20_prior[i] <- R_prior_i[idx_T20, idx_Res20]
rho150_prior[i] <- R_prior_i[idx_T150, idx_Res150]
}
cat("Resumen de correlaciones PRIOR simuladas:\n")## Resumen de correlaciones PRIOR simuladas:cat("rho(T20, Res20) — Media prior:", round(mean(rho20_prior), 4),
" SD:", round(sd(rho20_prior), 4), "\n")## rho(T20, Res20) — Media prior: -0.2967 SD: 0.1407cat("rho(T150, Res150) — Media prior:", round(mean(rho150_prior), 4),
" SD:", round(sd(rho150_prior), 4), "\n")## rho(T150, Res150) — Media prior: -0.2966 SD: 0.1411# -----------------------------------------------------------------------------
# 2. CONSTRUIR DATAFRAME PARA GGPLOT
# -----------------------------------------------------------------------------
# --- rho(T20, Res20) ---
df_comp_rho20 <- rbind(
data.frame(
rho = rho20_prior,
Distribucion = "Prior"
),
data.frame(
rho = rho20_post,
Distribucion = "Posterior"
)
)
# --- rho(T150, Res150) ---
df_comp_rho150 <- rbind(
data.frame(
rho = rho150_prior,
Distribucion = "Prior"
),
data.frame(
rho = rho150_post,
Distribucion = "Posterior"
)
)
# Convertir a factor con orden definido
df_comp_rho20$Distribucion <- factor(
df_comp_rho20$Distribucion,
levels = c("Prior", "Posterior")
)
df_comp_rho150$Distribucion <- factor(
df_comp_rho150$Distribucion,
levels = c("Prior", "Posterior")
)
# -----------------------------------------------------------------------------
# 3. LÍNEAS VERTICALES — medias de cada distribución
# -----------------------------------------------------------------------------
lineas_rho20 <- data.frame(
Distribucion = factor(c("Prior", "Posterior"),
levels = c("Prior", "Posterior")),
Media = c(mean(rho20_prior), mean(rho20_post))
)
lineas_rho150 <- data.frame(
Distribucion = factor(c("Prior", "Posterior"),
levels = c("Prior", "Posterior")),
Media = c(mean(rho150_prior), mean(rho150_post))
)
# -----------------------------------------------------------------------------
# 4. PALETA DE COLORES
# -----------------------------------------------------------------------------
colores <- c(
"Prior" = "#B0BEC5", # gris azulado — incertidumbre inicial
"Posterior" = "#2E86AB" # azul — actualización con datos
)
colores_linea <- c(
"Prior" = "#546E7A",
"Posterior" = "#E63946"
)
# -----------------------------------------------------------------------------
# 5. GRÁFICO rho(T20, Res20)
# -----------------------------------------------------------------------------
p_comp_rho20 <- ggplot(
df_comp_rho20,
aes(x = rho, fill = Distribucion, color = Distribucion)
) +
geom_density(
alpha = 0.45,
linewidth = 0.9
) +
geom_vline(
data = lineas_rho20,
aes(xintercept = Media, color = Distribucion),
linewidth = 1.1,
linetype = "dashed"
) +
scale_fill_manual(values = colores) +
scale_color_manual(values = colores_linea) +
annotate(
"text",
x = mean(rho20_prior) + 0.03,
y = Inf,
label = paste0("Media prior\n",
round(mean(rho20_prior), 3)),
vjust = 1.5,
hjust = 0,
size = 3.5,
color = colores_linea["Prior"]
) +
annotate(
"text",
x = mean(rho20_post) - 0.03,
y = Inf,
label = paste0("Media posterior\n",
round(mean(rho20_post), 3)),
vjust = 1.5,
hjust = 1,
size = 3.5,
color = colores_linea["Posterior"]
) +
labs(
title = expression(paste(
"Prior vs Posterior — ",
rho, "(T"[20], ", Res"[20], ")"
)),
subtitle = paste0(
"La distribución se concentra alrededor de ",
round(mean(rho20_post), 3),
" tras incorporar ", nrow(datos_std), " observaciones"
),
x = expression(rho),
y = "Densidad",
fill = "Distribución",
color = "Distribución"
) +
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(face = "bold", size = 13),
plot.subtitle = element_text(size = 10, color = "gray40"),
legend.position = "top"
)
# -----------------------------------------------------------------------------
# 6. GRÁFICO rho(T150, Res150)
# -----------------------------------------------------------------------------
p_comp_rho150 <- ggplot(
df_comp_rho150,
aes(x = rho, fill = Distribucion, color = Distribucion)
) +
geom_density(
alpha = 0.45,
linewidth = 0.9
) +
geom_vline(
data = lineas_rho150,
aes(xintercept = Media, color = Distribucion),
linewidth = 1.1,
linetype = "dashed"
) +
scale_fill_manual(values = colores) +
scale_color_manual(values = colores_linea) +
annotate(
"text",
x = mean(rho150_prior) + 0.03,
y = Inf,
label = paste0("Media prior\n",
round(mean(rho150_prior), 3)),
vjust = 1.5,
hjust = 0,
size = 3.5,
color = colores_linea["Prior"]
) +
annotate(
"text",
x = mean(rho150_post) - 0.03,
y = Inf,
label = paste0("Media posterior\n",
round(mean(rho150_post), 3)),
vjust = 1.5,
hjust = 1,
size = 3.5,
color = colores_linea["Posterior"]
) +
labs(
title = expression(paste(
"Prior vs Posterior — ",
rho, "(T"[150], ", Res"[150], ")"
)),
subtitle = paste0(
"La distribución se concentra alrededor de ",
round(mean(rho150_post), 3),
" tras incorporar ", nrow(datos_std), " observaciones"
),
x = expression(rho),
y = "Densidad",
fill = "Distribución",
color = "Distribución"
) +
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(face = "bold", size = 13),
plot.subtitle = element_text(size = 10, color = "gray40"),
legend.position = "top"
)
# -----------------------------------------------------------------------------
# 7. COMBINAR Y GUARDAR
# -----------------------------------------------------------------------------
fig_prior_posterior <- grid.arrange(
p_comp_rho20,
p_comp_rho150,
ncol = 1,
top = grid::textGrob(
"Actualización Bayesiana: Prior → Posterior",
gp = grid::gpar(fontface = "bold", fontsize = 14)
)
)
ggsave(
"fig_prior_vs_posterior.png",
fig_prior_posterior,
width = 10,
height = 10,
dpi = 200,
bg = "white"
)
cat("Figura guardada: fig_prior_vs_posterior.png\n")## Figura guardada: fig_prior_vs_posterior.png