# 1. Cargar librerías necesarias
library(lmtest)## Warning: package 'lmtest' was built under R version 4.5.3## Cargando paquete requerido: zoo##
## Adjuntando el paquete: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(sandwich)## Warning: package 'sandwich' was built under R version 4.5.3library(forecast)
# Definir la semilla global requerida
set.seed(50)# 1. Coeficiente de Determinación (R2)
R2 <- function(predicho, real) {
1 - (sum((real - predicho)^2) / sum((real - mean(real))^2))
}
# 2. Raíz del Error Cuadrático Medio (RMSE)
RMSE <- function(predicho, real) {
sqrt(mean((real - predicho)^2))
}
# 3. Error Absoluto Medio (MAE)
MAE <- function(predicho, real) {
mean(abs(real - predicho))
}
# 4. Error Porcentual Absoluto Medio (MAPE)
MAPE <- function(predicho, real) {
mean(abs((real - predicho) / real))
}
# 5. U de Theil (Tipo 1)
TheilU <- function(predicho, real, type = 1) {
sqrt(mean((predicho - real)^2)) / (sqrt(mean(predicho^2)) + sqrt(mean(real^2)))
}
# 6. Proporción de Sesgo (Um)
Um <- function(predicho, real) {
(mean(predicho) - mean(real))^2 / mean((predicho - real)^2)
}
# 7. Proporción de Varianza (Us)
Us <- function(predicho, real) {
(sd(predicho) - sd(real))^2 / mean((predicho - real)^2)
}
# 8. Proporción de Covarianza (Uc)
Uc <- function(predicho, real) {
2 * (1 - cor(predicho, real)) * sd(predicho) * sd(real) / mean((predicho - real)^2)
}#configuracion inicial de la simulación, incluyendo la semilla global para garantizar la reproducibilidad de los resultados.
library(mlbench)## Warning: package 'mlbench' was built under R version 4.5.3library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.5.3## Warning: package 'dplyr' was built under R version 4.5.3## Warning: package 'lubridate' was built under R version 4.5.3## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.1 ✔ readr 2.2.0
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.2 ✔ tibble 3.3.1
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
## ✔ purrr 1.2.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsset.seed(50)
data(BostonHousing)
numero_de_muestras <- 100 # Cantidad de réplicas de la presentación
proporcion_entrenamiento <- 0.75 # 75% para entrenamientomuestras <- list()
for(i in 1:numero_de_muestras) {
muestras[[i]] <- sample(1:nrow(BostonHousing), size = floor(proporcion_entrenamiento * nrow(BostonHousing)))
}Modelos_Entrenamiento <- vector(mode = "list", length = numero_de_muestras)
Pronostico_Prueba <- vector(mode = "list", length = numero_de_muestras)
Resultados_Performance_data_entrenamiento <- vector(mode = "list", length = numero_de_muestras)
Resultados_Performance <- vector(mode = "list", length = numero_de_muestras)for(j in 1:numero_de_muestras){
# División de datos según los índices guardados en 'muestras'
Datos_Entrenamiento <- BostonHousing[muestras[[j]], ]
Datos_Prueba <- BostonHousing[-muestras[[j]], ]
# Estimación del modelo lineal para la muestra actual
Modelos_Entrenamiento[[j]] <- lm(formula = medv ~ ., data = Datos_Entrenamiento)
# Pronóstico con los datos de prueba (datos "desconocidos")
Pronostico_Prueba[[j]] <- Modelos_Entrenamiento[[j]] %>% predict(Datos_Prueba)
Resultados_Performance_data_entrenamiento[[j]] <- data.frame(
R2 = R2(Modelos_Entrenamiento[[j]]$fitted.values, Datos_Entrenamiento$medv),
RMSE = RMSE(Modelos_Entrenamiento[[j]]$fitted.values, Datos_Entrenamiento$medv),
MAE = MAE(Modelos_Entrenamiento[[j]]$fitted.values, Datos_Entrenamiento$medv),
MAPE = MAPE(Modelos_Entrenamiento[[j]]$fitted.values, Datos_Entrenamiento$medv) * 100,
THEIL = TheilU(Modelos_Entrenamiento[[j]]$fitted.values, Datos_Entrenamiento$medv, type = 1),
Um = Um(Modelos_Entrenamiento[[j]]$fitted.values, Datos_Entrenamiento$medv),
Us = Us(Modelos_Entrenamiento[[j]]$fitted.values, Datos_Entrenamiento$medv),
Uc = Uc(Modelos_Entrenamiento[[j]]$fitted.values, Datos_Entrenamiento$medv)
)
# Cálculo de métricas en Datos de Prueba (Pronóstico)
Resultados_Performance[[j]] <- data.frame(
R2 = cor(Pronostico_Prueba[[j]], Datos_Prueba$medv)^2, # R² elemental
RMSE = sqrt(mean((Datos_Prueba$medv - Pronostico_Prueba[[j]])^2)),
MAE = mean(abs(Datos_Prueba$medv - Pronostico_Prueba[[j]])),
MAPE = mean(abs((Datos_Prueba$medv - Pronostico_Prueba[[j]]) / Datos_Prueba$medv)) * 100
)
} # Fin del bucle # 6. Consolidación de todos los resultados en DataFrames
df_performance_entrenamiento <- do.call(rbind, Resultados_Performance_data_entrenamiento)
df_performance_prueba <- do.call(rbind, Resultados_Performance)# 7. Visualización del resumen estadístico de la simulación (Datos de Prueba)
print("--- RESUMEN DE RENDIMIENTO EN DATOS DE PRUEBA ---")## [1] "--- RESUMEN DE RENDIMIENTO EN DATOS DE PRUEBA ---"summary(df_performance_prueba)## R2 RMSE MAE MAPE
## Min. :0.5587 Min. :3.377 Min. :2.679 Min. :13.69
## 1st Qu.:0.6906 1st Qu.:4.465 1st Qu.:3.201 1st Qu.:15.94
## Median :0.7289 Median :4.768 Median :3.403 Median :16.99
## Mean :0.7253 Mean :4.823 Mean :3.392 Mean :17.06
## 3rd Qu.:0.7604 3rd Qu.:5.198 3rd Qu.:3.576 3rd Qu.:18.22
## Max. :0.8317 Max. :6.239 Max. :4.004 Max. :21.40#Cargar librería y traer los datos del paquete de origen
library(tidyverse)
library(sandwich)
library(wooldridge) ## Warning: package 'wooldridge' was built under R version 4.5.3data("smoke")
datos_smoke <- smoke
n_total <- nrow(datos_smoke) # 2. Configuración de parámetros de la simulación
set.seed(50)
replicas <- 500
proporcion_train <- 0.75
n_total <- nrow(datos_smoke) # 3. Generar las submuestras aleatorias para las 500 réplicas
muestras_hac <- list()
for(i in 1:replicas) {
muestras_hac[[i]] <- sample(1:n_total, size = floor(proporcion_train * n_total))
}# 4. Inicializar listas para almacenar los resultados
Modelos_Entrenamiento_HAC <- vector(mode = "list", length = replicas)
Pronostico_Prueba_HAC <- vector(mode = "list", length = replicas)
Resultados_Performance_HAC <- vector(mode = "list", length = replicas)# 5. Bucle de simulación para las 500 réplicas
for(j in 1:replicas) {
# División de conjuntos (75% Entrenamiento / 25% Prueba)
Datos_Entrenamiento <- datos_smoke[muestras_hac[[j]], ]
Datos_Prueba <- datos_smoke[-muestras_hac[[j]], ]
# Estimación del modelo especificado en tu tarea anterior
Modelos_Entrenamiento_HAC[[j]] <- lm(
cigs ~ cigpric + lcigpric + income + lincome + age + agesq + educ + white + restaurn,
data = Datos_Entrenamiento
)
# Pronóstico sobre los datos "desconocidos" (Datos_Prueba)
Pronostico_Prueba_HAC[[j]] <- Modelos_Entrenamiento_HAC[[j]] %>% predict(Datos_Prueba)
# Cálculo de métricas de desempeño con los datos de prueba
Resultados_Performance_HAC[[j]] <- data.frame(
R2 = R2(Pronostico_Prueba_HAC[[j]], Datos_Prueba$cigs),
RMSE = RMSE(Pronostico_Prueba_HAC[[j]], Datos_Prueba$cigs),
MAE = MAE(Pronostico_Prueba_HAC[[j]], Datos_Prueba$cigs),
MAPE = MAPE(Pronostico_Prueba_HAC[[j]], Datos_Prueba$cigs) * 100,
THEIL = TheilU(Pronostico_Prueba_HAC[[j]], Datos_Prueba$cigs, type = 1),
Um = Um(Pronostico_Prueba_HAC[[j]], Datos_Prueba$cigs),
Us = Us(Pronostico_Prueba_HAC[[j]], Datos_Prueba$cigs),
Uc = Uc(Pronostico_Prueba_HAC[[j]], Datos_Prueba$cigs)
)
}# 6. Consolidación de todos los resultados en un DataFrame
df_performance_hac <- do.call(rbind, Resultados_Performance_HAC)#Analisis de los resultados del desempeño del modelo
#Obtención del resumen estadístico general
print("--- RESUMEN ESTADÍSTICO DE DESEMPEÑO (500 RÉPLICAS) ---")## [1] "--- RESUMEN ESTADÍSTICO DE DESEMPEÑO (500 RÉPLICAS) ---"summary(df_performance_hac)## R2 RMSE MAE MAPE
## Min. :-0.15478 Min. :11.20 Min. : 9.608 Min. :Inf
## 1st Qu.: 0.01170 1st Qu.:12.86 1st Qu.:10.412 1st Qu.:Inf
## Median : 0.03080 Median :13.48 Median :10.714 Median :Inf
## Mean : 0.02794 Mean :13.49 Mean :10.723 Mean :Inf
## 3rd Qu.: 0.04675 3rd Qu.:14.11 3rd Qu.:11.029 3rd Qu.:Inf
## Max. : 0.09272 Max. :16.00 Max. :12.031 Max. :Inf
## THEIL Um Us Uc
## Min. :0.4859 Min. :1.320e-07 Min. :0.3287 Min. :0.2352
## 1st Qu.:0.5168 1st Qu.:6.507e-04 1st Qu.:0.5592 1st Qu.:0.3634
## Median :0.5279 Median :2.835e-03 Median :0.6034 Median :0.3983
## Mean :0.5282 Mean :6.842e-03 Mean :0.5960 Mean :0.4021
## 3rd Qu.:0.5400 3rd Qu.:8.419e-03 3rd Qu.:0.6372 3rd Qu.:0.4391
## Max. :0.5745 Max. :1.019e-01 Max. :0.7437 Max. :0.6013#Promedio de las proporciones de error de la U de Theil
print("--- PROMEDIOS DE DESCOMPOSICIÓN DE THEIL ---")## [1] "--- PROMEDIOS DE DESCOMPOSICIÓN DE THEIL ---"colMeans(df_performance_hac[, c("Um", "Us", "Uc")])## Um Us Uc
## 0.006842415 0.596007411 0.402091257