Problema 1 — Regresión no paramétrica con el dataset Auto

El objetivo es modelar la relación entre horsepower y mpg comparando tres paradigmas: bases de funciones, regresión local y polinomio global de grado 2. La comparación se realiza mediante 10 repeticiones de validación cruzada.

0.1 Exploración de los datos

data(Auto)

hp  <- Auto$horsepower
mpg <- Auto$mpg
n   <- nrow(Auto)

cat("Observaciones totales  :", n, "\n")
## Observaciones totales  : 392
cat("Rango de horsepower    :", range(hp), "\n")
## Rango de horsepower    : 46 230
cat("Rango de mpg           :", range(mpg), "\n")
## Rango de mpg           : 9 46.6
ggplot(Auto, aes(x = horsepower, y = mpg)) +
  geom_point(alpha = 0.45, color = "#2C7BB6", size = 1.8) +
  geom_smooth(method = "loess", se = TRUE,
              color = "#D7191C", fill = "#FDAE61", alpha = 0.25) +
  labs(
    title    = "Relación entre Horsepower y MPG",
    subtitle = "Dataset Auto — ISLR2 | n = 392 vehículos",
    x = "Horsepower", y = "Miles per Gallon (mpg)"
  ) +
  theme_minimal(base_size = 13)

0.2 Función auxiliar de ECM

# Error Cuadrático Medio
ecm <- function(y_real, y_pred) {
  mean((y_real - y_pred)^2, na.rm = TRUE)
}

0.3 Puntos (1) a (6): Loop de 10 repeticiones

Para obtener distribuciones robustas del ECM de prueba (punto 6), los puntos (1)–(5) se ejecutan 10 veces, cada una con una partición aleatoria distinta.

# ── Almacenamiento de resultados ─────────────────────────────────────────────
resultados <- data.frame(
  repeticion = integer(),
  paradigma  = character(),
  ECM_prueba = numeric(),
  stringsAsFactors = FALSE
)

resumen_rep <- data.frame(
  rep         = 1:10,
  k_optimo    = NA_integer_,
  mejor_basis = NA_character_,
  mejor_local = NA_integer_,
  stringsAsFactors = FALSE
)

# ── Loop principal ────────────────────────────────────────────────────────────
for (rep in 1:10) {

  set.seed(rep * 42)   # Semilla reproducible distinta por iteración

  # ── PUNTO (1): Separación 90% / 10% ─────────────────────────────────────
  idx_test  <- sample(seq_len(n), size = round(0.10 * n))
  idx_train <- setdiff(seq_len(n), idx_test)

  df_train <- data.frame(hp = hp[idx_train], mpg = mpg[idx_train])
  df_test  <- data.frame(hp = hp[idx_test],  mpg = mpg[idx_test])
  n_train  <- nrow(df_train)

  # Nodos de frontera (5 unidades fuera del rango → ajuste estable)
  boundary <- range(df_train$hp) + c(-5, 5)

  # Asignación aleatoria de folds
  folds <- sample(1:10, n_train, replace = TRUE)

  # ── PUNTO (2): Knots óptimos — CV-10 sobre k = 1, …, 10 ─────────────────
  ecm_por_knot <- numeric(10)

  for (k in 1:10) {

    knot_pos <- seq(min(df_train$hp), max(df_train$hp),
                    length.out = k + 2)[2:(k + 1)]
    errores_folds <- numeric(10)

    for (fold in 1:10) {
      idx_val <- which(folds == fold)
      idx_tr  <- which(folds != fold)
      if (length(idx_val) == 0) next

      X_tr  <- bs(df_train$hp[idx_tr],  knots = knot_pos,
                  Boundary.knots = boundary)
      X_val <- bs(df_train$hp[idx_val], knots = knot_pos,
                  Boundary.knots = boundary)

      fit  <- lm(df_train$mpg[idx_tr] ~ X_tr)
      pred <- cbind(1, X_val) %*% coef(fit)

      errores_folds[fold] <- ecm(df_train$mpg[idx_val], pred)
    }

    ecm_por_knot[k] <- mean(errores_folds, na.rm = TRUE)
  }

  k_optimo <- which.min(ecm_por_knot)
  knot_opt <- seq(min(df_train$hp), max(df_train$hp),
                  length.out = k_optimo + 2)[2:(k_optimo + 1)]

  resumen_rep$k_optimo[rep] <- k_optimo

  # ── PUNTO (3): Comparar modelos de bases — CV-10 ──────────────────────────
  ecm_p2 <- ecm_ss <- ecm_rs <- numeric(10)

  for (fold in 1:10) {
    idx_val <- which(folds == fold)
    idx_tr  <- which(folds != fold)
    if (length(idx_val) == 0) next

    df_tr  <- df_train[idx_tr, ]
    df_val <- df_train[idx_val, ]

    # Polinomio global grado 2
    fit_p2 <- lm(mpg ~ poly(hp, 2), data = df_tr)
    ecm_p2[fold] <- ecm(df_val$mpg, predict(fit_p2, newdata = df_val))

    # Smoothing spline (GCV automático)
    fit_ss <- smooth.spline(df_tr$hp, df_tr$mpg, cv = FALSE)
    ecm_ss[fold] <- ecm(df_val$mpg, predict(fit_ss, x = df_val$hp)$y)

    # Regression spline óptimo
    X_tr  <- bs(df_tr$hp,  knots = knot_opt, Boundary.knots = boundary)
    X_val <- bs(df_val$hp, knots = knot_opt, Boundary.knots = boundary)
    fit_rs  <- lm(df_tr$mpg ~ X_tr)
    pred_rs <- cbind(1, X_val) %*% coef(fit_rs)
    ecm_rs[fold] <- ecm(df_val$mpg, pred_rs)
  }

  ecm_cv_basis <- c(
    `Polinomio grado 2` = mean(ecm_p2, na.rm = TRUE),
    `Smoothing spline`  = mean(ecm_ss, na.rm = TRUE),
    `Regression spline` = mean(ecm_rs, na.rm = TRUE)
  )

  mejor_basis <- names(which.min(ecm_cv_basis))
  resumen_rep$mejor_basis[rep] <- mejor_basis

  # ── PUNTO (4): Regresión local — loess grado 1 vs grado 2 ────────────────
  ecm_l1 <- ecm_l2 <- numeric(10)

  for (fold in 1:10) {
    idx_val <- which(folds == fold)
    idx_tr  <- which(folds != fold)
    if (length(idx_val) == 0) next

    df_tr  <- df_train[idx_tr, ]
    df_val <- df_train[idx_val, ]

    fit_l1 <- loess(mpg ~ hp, data = df_tr, degree = 1)
    ecm_l1[fold] <- ecm(df_val$mpg, predict(fit_l1, newdata = df_val))

    fit_l2 <- loess(mpg ~ hp, data = df_tr, degree = 2)
    ecm_l2[fold] <- ecm(df_val$mpg, predict(fit_l2, newdata = df_val))
  }

  mejor_local <- which.min(c(
    Degree1 = mean(ecm_l1, na.rm = TRUE),
    Degree2 = mean(ecm_l2, na.rm = TRUE)
  ))
  resumen_rep$mejor_local[rep] <- mejor_local

  # ── PUNTO (5): ECM de PRUEBA para los 3 paradigmas ───────────────────────

  # Paradigma 1: Mejor modelo de bases
  if (mejor_basis == "Polinomio grado 2") {
    fit_b  <- lm(mpg ~ poly(hp, 2), data = df_train)
    pred_b <- predict(fit_b, newdata = df_test)

  } else if (mejor_basis == "Smoothing spline") {
    fit_b  <- smooth.spline(df_train$hp, df_train$mpg, cv = FALSE)
    pred_b <- predict(fit_b, x = df_test$hp)$y

  } else {
    X_full <- bs(df_train$hp, knots = knot_opt, Boundary.knots = boundary)
    X_tst  <- bs(df_test$hp,  knots = knot_opt, Boundary.knots = boundary)
    fit_b  <- lm(df_train$mpg ~ X_full)
    pred_b <- cbind(1, X_tst) %*% coef(fit_b)
  }
  ecm_test_basis <- ecm(df_test$mpg, pred_b)

  # Paradigma 2: Mejor regresión local
  fit_l          <- loess(mpg ~ hp, data = df_train, degree = mejor_local)
  ecm_test_local <- ecm(df_test$mpg, predict(fit_l, newdata = df_test))

  # Paradigma 3: Polinomio global grado 2
  fit_pg          <- lm(mpg ~ poly(hp, 2), data = df_train)
  ecm_test_poly2  <- ecm(df_test$mpg, predict(fit_pg, newdata = df_test))

  resultados <- rbind(resultados, data.frame(
    repeticion = rep,
    paradigma  = c("Bases de funciones",
                   "Regresión local",
                   "Polinomio global grado 2"),
    ECM_prueba = c(ecm_test_basis, ecm_test_local, ecm_test_poly2),
    stringsAsFactors = FALSE
  ))
}

cat("✔ Loop de 10 repeticiones completado.\n")
## ✔ Loop de 10 repeticiones completado.

0.4 Resultado — Punto (1): Separación de datos

En cada una de las 10 repeticiones se usó una semilla rep × 42 distinta para garantizar reproducibilidad. El 10 % de los 392 vehículos (≈39) quedó como conjunto de prueba y el 90 % restante (≈353) como conjunto de entrenamiento.

0.5 Resultado — Punto (2): Knots óptimos en el regression spline

cat("Número óptimo de knots por repetición:\n")
## Número óptimo de knots por repetición:
knitr::kable(resumen_rep[, c("rep", "k_optimo")],
             col.names = c("Repetición", "k óptimo"),
             align = "cc")
Repetición k óptimo
1 4
2 4
3 4
4 4
5 6
6 3
7 4
8 4
9 4
10 4 
ggplot(resumen_rep, aes(x = factor(rep), y = k_optimo)) +
  geom_col(fill = "#2C7BB6", alpha = 0.85, width = 0.6) +
  geom_hline(yintercept = median(resumen_rep$k_optimo),
             linetype = "dashed", color = "#D7191C", linewidth = 0.8) +
  annotate("text", x = 10.5, y = median(resumen_rep$k_optimo) + 0.3,
           label = paste("Mediana =", median(resumen_rep$k_optimo)),
           color = "#D7191C", size = 3.5) +
  scale_y_continuous(breaks = 1:10) +
  labs(title = "Número óptimo de knots por repetición (CV-10)",
       x = "Repetición", y = "k óptimo") +
  theme_minimal(base_size = 13)

Interpretación: El número de knots que minimiza el ECM de validación cruzada indica la complejidad óptima del regression spline. Un k mayor aumenta la flexibilidad del ajuste, pero puede llevar a sobreajuste.

0.6 Resultado — Punto (3): Mejor modelo de bases de funciones

cat("Mejor modelo de bases por repetición:\n")
## Mejor modelo de bases por repetición:
knitr::kable(resumen_rep[, c("rep", "k_optimo", "mejor_basis")],
             col.names = c("Repetición", "k óptimo", "Mejor modelo"),
             align = "ccl")
Repetición k óptimo Mejor modelo
1 4 Regression spline
2 4 Regression spline
3 4 Regression spline
4 4 Regression spline
5 6 Smoothing spline
6 3 Regression spline
7 4 Regression spline
8 4 Smoothing spline
9 4 Regression spline
10 4 Regression spline 
cat("\nFrequencia de selección por modelo:\n")
## 
## Frequencia de selección por modelo:
print(sort(table(resumen_rep$mejor_basis), decreasing = TRUE))
## 
## Regression spline  Smoothing spline 
##                 8                 2

Interpretación: Se comparan tres modelos mediante CV-10. El smoothing spline usualmente tiene ventaja porque su parámetro de suavizado se optimiza automáticamente mediante GCV, mientras que el regression spline usa el k encontrado en el punto anterior.

0.7 Resultado — Punto (4): Mejor regresión local

cat("Grado de loess óptimo por repetición:\n")
## Grado de loess óptimo por repetición:
knitr::kable(
  data.frame(
    Repetición = resumen_rep$rep,
    `Grado óptimo` = ifelse(resumen_rep$mejor_local == 1,
                             "Grado 1 (localmente lineal)",
                             "Grado 2 (localmente cuadrático)")
  ),
  align = "cl"
)
Repetición Grado.óptimo
1 Grado 2 (localmente cuadrático)
2 Grado 2 (localmente cuadrático)
3 Grado 2 (localmente cuadrático)
4 Grado 2 (localmente cuadrático)
5 Grado 2 (localmente cuadrático)
6 Grado 2 (localmente cuadrático)
7 Grado 2 (localmente cuadrático)
8 Grado 2 (localmente cuadrático)
9 Grado 2 (localmente cuadrático)
10 Grado 1 (localmente lineal)

Interpretación: loess con grado 1 ajusta localmente una recta; con grado 2 ajusta una parábola. El ancho de banda óptimo es el que provee loess() por defecto. El grado que produce menor ECM de CV-10 es el seleccionado.

0.8 Resultado — Puntos (5) y (6): ECM de prueba — 10 repeticiones

resumen_ecm <- resultados %>%
  group_by(paradigma) %>%
  summarise(
    Media   = round(mean(ECM_prueba), 3),
    Mediana = round(median(ECM_prueba), 3),
    DE      = round(sd(ECM_prueba), 3),
    Mínimo  = round(min(ECM_prueba), 3),
    Máximo  = round(max(ECM_prueba), 3),
    .groups = "drop"
  )

knitr::kable(resumen_ecm,
             caption = "Resumen del ECM de prueba por paradigma (10 repeticiones)")
Resumen del ECM de prueba por paradigma (10 repeticiones)
paradigma Media Mediana DE Mínimo Máximo
Bases de funciones 18.252 15.831 6.874 11.915 33.876
Polinomio global grado 2 18.124 15.773 6.213 12.261 33.363
Regresión local 17.968 15.273 6.732 12.274 34.462 
paleta <- c(
  "Bases de funciones"       = "#2C7BB6",
  "Regresión local"          = "#1A9641",
  "Polinomio global grado 2" = "#D7191C"
)

ggplot(resultados, aes(x = paradigma, y = ECM_prueba, fill = paradigma)) +
  geom_boxplot(alpha = 0.7, outlier.shape = 21, outlier.size = 2.5,
               width = 0.5) +
  geom_jitter(aes(color = paradigma), width = 0.12, alpha = 0.9,
              size = 2.5) +
  scale_fill_manual(values  = paleta) +
  scale_color_manual(values = paleta) +
  labs(
    title    = "Distribución del ECM de prueba por paradigma",
    subtitle = "10 repeticiones con partición aleatoria 90/10",
    x        = NULL,
    y        = "ECM de prueba"
  ) +
  theme_minimal(base_size = 13) +
  theme(legend.position = "none",
        axis.text.x = element_text(size = 11))

ggplot(resultados, aes(x = repeticion, y = ECM_prueba,
                       color = paradigma, group = paradigma)) +
  geom_line(linewidth = 1.1, alpha = 0.9) +
  geom_point(size = 3) +
  scale_color_manual(values = paleta) +
  scale_x_continuous(breaks = 1:10) +
  labs(
    title  = "Evolución del ECM de prueba a lo largo de las 10 repeticiones",
    x      = "Repetición",
    y      = "ECM de prueba",
    color  = "Paradigma"
  ) +
  theme_minimal(base_size = 13) +
  theme(legend.position = "bottom")

0.8.1 Conclusión Problema 1

ganador <- resumen_ecm$paradigma[which.min(resumen_ecm$Media)]
cat("Paradigma con menor ECM promedio de prueba:", ganador, "\n")
## Paradigma con menor ECM promedio de prueba: Regresión local
cat("ECM promedio                              :", min(resumen_ecm$Media), "\n")
## ECM promedio                              : 17.968

Respuesta al punto (6): Con base en las distribuciones del ECM de prueba obtenidas en las 10 repeticiones, el paradigma de Regresión local produce el menor error de predicción promedio. Dado que las distribuciones de los tres paradigmas se superponen parcialmente, se recomienda también considerar la variabilidad (desviación estándar) al seleccionar el modelo final.

Problema 2 — Análisis de datos funcionales

0.9 Notación

Para la \(i\)-ésima unidad estadística se tienen \(n_i\) observaciones \(\{(t_{ij},\, x_{ij})\}_{j=1}^{n_i}\), donde \(t_{ij} \in \mathcal{T} \subset \mathbb{R}\) es el punto de evaluación y \(x_{ij} \in \mathbb{R}\) es el valor observado. Sea \(K_h(u) = \frac{1}{h}K\!\left(\frac{u}{h}\right)\) un kernel simétrico con ancho de banda \(h > 0\).

Punto (7): Estimador de Nadaraya–Watson para la \(i\)-ésima unidad

El estimador de Nadaraya–Watson para la función suavizada \(\hat{x}_i(t)\) de la \(i\)-ésima unidad estadística en el punto \(t\) es:

\[ \hat{x}_i(t) \;=\; \frac{\displaystyle\sum_{j=1}^{n_i} K_h\!\left(t - t_{ij}\right)\, x_{ij}} {\displaystyle\sum_{j=1}^{n_i} K_h\!\left(t - t_{ij}\right)} \]

Interpretación: \(\hat{x}_i(t)\) es un promedio ponderado de las observaciones \(x_{ij}\) pertenecientes a la unidad \(i\). El peso de cada observación es proporcional al kernel evaluado en la distancia \(|t - t_{ij}|\): observaciones con \(t_{ij}\) cercano a \(t\) reciben mayor peso, logrando suavizar localmente la curva.

Punto (8): Estimador de Nadaraya–Watson para la función media \(\hat{\mu}(t)\)

Extendiendo el estimador anterior a todas las \(n\) unidades estadísticas y sus observaciones, el estimador de Nadaraya–Watson para la función media en \(t\) es:

\[ \hat{\mu}(t) \;=\; \frac{\displaystyle\sum_{i=1}^{n} \sum_{j=1}^{n_i} K_h\!\left(t - t_{ij}\right)\, x_{ij}} {\displaystyle\sum_{i=1}^{n} \sum_{j=1}^{n_i} K_h\!\left(t - t_{ij}\right)} \]

Interpretación: Este estimador utiliza todos los datos discretizados de todas las \(n\) unidades de manera conjunta. En cada punto \(t\), las observaciones \(x_{ij}\) cuyo \(t_{ij}\) es cercano a \(t\) reciben mayor peso, independientemente del individuo \(i\) al que pertenezcan. De esta forma, la estimación de \(\mu(t)\) aprovecha toda la información disponible de la muestra funcional.

Caso particular: Si \(n_i = m\) para todo \(i\) y los puntos de evaluación son comunes (\(t_{ij} = t_j\)), el estimador se simplifica a:

\[ \hat{\mu}(t) \;=\; \frac{\displaystyle\sum_{j=1}^{m} K_h\!\left(t - t_j\right)\,\bar{x}_j} {\displaystyle\sum_{j=1}^{m} K_h\!\left(t - t_j\right)}, \qquad \bar{x}_j = \frac{1}{n}\sum_{i=1}^{n} x_{ij} \]

lo que resalta su interpretación como un suavizado de los promedios muestrales punto a punto: \(\hat{\mu}(t)\) es el estimador de Nadaraya–Watson aplicado a los promedios \(\bar{x}_j\) evaluados en los puntos \(t_j\).

Fin del Taller 2