Modelos Estadísticos I

Regresión no paramétrica

library(dplyr)       
library(ggplot2)
library(MASS)
library(mgcv)
library(nortest)

1 Modelos paramétricos

1.1 Modelo lineal simple

En el caso usual, el paramétrico, la regresión lineal simple según Ochoa (2020a) es uno de los modelos más comunes. Un modelo estadístico lineal general se escribe como:

\[ y = f(x) + \varepsilon \tag{1}\label{eq:modelo} \]

donde \(y\) y \(\varepsilon\) son variables aleatorias, y \(f\) es una función determinística de una variable no aleatoria.
En un modelo lineal, \(f\) es lineal en los parámetros desconocidos:

\[ f(x) = \beta_0 + \beta_1 x \]

Así, el modelo se expresa como:

\[ y = \beta_0 + \beta_1 x + \varepsilon, \quad \mathbb{E}[\varepsilon] = 0, \quad \text{Var}(\varepsilon) = \sigma^2 \]

Los estimadores de mínimos cuadrados para \(\beta_0\) y \(\beta_1\) son:

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

El modelo estimado es entonces:

\[ \hat{f}(x) = \hat{\beta}_0 + \hat{\beta}_1 x \]

Para los siguientes datos simulados, se calcula y grafica la regresión de los datos:

set.seed(123)
x <- seq(1, 10, length.out = 100)
y <- exp(0.3 * x) + rgamma(100, shape = 2, rate = 0.5)  
# Ajuste lineal
model <- lm(y ~ x)
summary(model)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.6518 -2.0117 -0.7275  1.8043  9.9492 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.6421     0.6961  -0.922    0.359    
## x             2.0113     0.1142  17.609   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.997 on 98 degrees of freedom
## Multiple R-squared:  0.7598, Adjusted R-squared:  0.7574 
## F-statistic: 310.1 on 1 and 98 DF,  p-value: < 2.2e-16

# Opción 1: Usar broom para un formato tidy
library(broom)
tidy_summary <- tidy(summary(model))
knitr::kable(tidy_summary, digits = 3)

Table 1.1: Summary del modelo ajustado

term	estimate	std.error	statistic	p.value
(Intercept)	-0.642	0.696	-0.922	0.359
x	2.011	0.114	17.609	0.000

library(ggplot2)
ggplot(data.frame(x, y), aes(x, y)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  ggtitle("Ajuste lineal")

## `geom_smooth()` using formula = 'y ~ x'

Figure 1.1: Ajuste lineal simple

Como se puede observar en 1.1 y en 1.1 el modelo para ser no muy adecuado, esto se discutira mas adelante.

1.2 Modelos lineal múltiple

Supongamos que tenemos datos \((x_1, Y_1), \ldots, (x_n, Y_n)\) donde \(Y_i \in \mathbb{R}\) y \(x_i = (x_{i1}, \ldots, x_{ip})^T \in \mathbb{R}^p\). El modelo de regresión lineal asume que

\[ Y_i = r(x_i) + \varepsilon_i \equiv \sum_{j=1}^{p} \beta_j x_{ij} + \varepsilon_i, \quad i = 1, \ldots, n, (*) \]

donde \(\mathbb{E}(\varepsilon_i) = 0\) y \(\mathbb{V}(\varepsilon_i) = \sigma^2\).

Usualmente, queremos incluir una ordenada al origen en el modelo, por lo que adoptaremos la convención de que \(x_{i1} = 1\).

La matriz de diseño \(X\) es la matriz \(n \times p\) definida por

\[ X = \begin{pmatrix} x_{11} & x_{12} & \cdots & x_{1p} \\ x_{21} & x_{22} & \cdots & x_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ x_{n1} & x_{n2} & \cdots & x_{np} \end{pmatrix}. \]

El conjunto \(\mathcal{L}\) de vectores que se pueden obtener como combinaciones lineales de las columnas de \(X\), se denomina el espacio columna de \(X\).

Sea \(Y = (Y_1, \ldots, Y_n)^T\), \(\varepsilon = (\varepsilon_1, \ldots, \varepsilon_n)^T\), y \(\beta = (\beta_1, \ldots, \beta_p)^T\). Podemos entonces escribir (*) como

\[ Y = X\beta + \varepsilon. \tag{5.5} \]

El estimador de mínimos cuadrados \(\hat{\beta} = (\hat{\beta}_1, \ldots, \hat{\beta}_p)^T\) es el vector que minimiza la suma de cuadrados de los residuos

\[ \text{RSS} = (Y - X\beta)^T (Y - X\beta) = \sum_{i=1}^{n} \left( Y_i - \sum_{j=1}^{p} x_{ij} \beta_j \right)^2. \]

Si el problema considera un predictor \(q\)-dimensional y existe \((\mathbf{X}^T\mathbf{X})^{-1}\), entonces el estimador de mínimos cuadrados \(\hat{\boldsymbol{\beta}}\) del vector de parámetros \(\boldsymbol{\beta}\) está dado por:

\[ \hat{\boldsymbol{\beta}} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{y} \]

donde \(\mathbf{X}\) y \(\mathbf{y}\) han sido definidos previamente para el modelo lineal general bajo el supuesto \(\mathbb{E}(\boldsymbol{\varepsilon}) = 0\), \(\hat{\boldsymbol{\beta}}\) es un estimador insesgado de \(\boldsymbol{\beta}\).

Por lo tanto para los regresores \(x\), \(x^2\), tenemos que:

set.seed(123)
x <- seq(1, 10, length.out = 100)
y <- exp(0.3 * x) + rgamma(100, shape = 2, rate = 0.5)  
# Ajuste lineal
model_cuadrado <- lm(y ~ poly(x,defree=2))
df <- data.frame(x = x, y = y)
ggplot(df, aes(x, y)) +
  geom_point() +
  geom_line(aes(y = predict(model_cuadrado)), color = "red") +
  ggtitle("Ajuste lineal con un polinomio  de grado 2")

summary(model_cuadrado)

## 
## Call:
## lm(formula = y ~ poly(x, defree = 2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.8361 -1.7778 -0.1533  1.0525  7.6576 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           10.4204     0.2228   46.78  < 2e-16 ***
## poly(x, defree = 2)1  52.7816     2.2275   23.70  < 2e-16 ***
## poly(x, defree = 2)2  19.9797     2.2275    8.97 2.25e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.228 on 97 degrees of freedom
## Multiple R-squared:  0.8687, Adjusted R-squared:  0.866 
## F-statistic:   321 on 2 and 97 DF,  p-value: < 2.2e-16

Así tenemos un modelo que para ser mas adecuado, no obstante:

par(mfrow=c(2,2),oma=c(0,0,0,0))
plot(model_cuadrado)

• Lo que se confirma con:

ad_test <- ad.test(model_cuadrado$residuals)
ad_test

## 
##  Anderson-Darling normality test
## 
## data:  model_cuadrado$residuals
## A = 1.8024, p-value = 0.0001216

1.3 Supuestos Clásicos de Regresión Lineal

La regresión lineal convencional se basa en dos pilares fundamentales:

1. Adecuación del modelo lineal: La relación verdadera es de naturaleza lineal
2. Propiedades de los errores:
   • Independencia: \(\varepsilon_i \perp \varepsilon_j\) para \(i \neq j\)
   • Distribución normal: \(\varepsilon_i \sim N(0, \sigma^2)\)
   • Homocedasticidad: \(Var(\varepsilon_i) = \sigma^2\) constante

1.4 Consecuencias de Violar los Supuestos

Cuando se ignoran estos supuestos, los procedimientos inferenciales tradicionales pierden validez, pues las propiedades óptimas de los estimadores, obtenidas mediante derivación de la función de verosimilitud \(L\), dependen críticamente de la normalidad

plot(model, which = c(1))

Figure 1.2: Supuestos en la regresión-1

plot(model, which = c(2))

Figure 1.3: Supuestos en la regresión-2

plot(model, which = c(3))

Figure 1.4: Supuestos en la regresión-3

ad_test <- ad.test(model$residuals)
ad_test

## 
##  Anderson-Darling normality test
## 
## data:  model$residuals
## A = 1.3389, p-value = 0.001704

Table 1.2: Resultados del test Anderson-Darling

	Estadístico_AD	Valor_p
A	1.338856	0.0017036

En base a la Figura 1.2,1.3,1.4, se observa que los supuestos del modelo no se cumplen adecuadamente, especialmente la normalidad de los residuos (como lo indica el test de Anderson-Darling) (Tabla 1.2) Por lo tanto, se opta por aplicar una transformación Box-Cox, con el objetivo de mejorar el ajuste del modelo y aproximarse al cumplimiento de los supuestos.

Así aplicando el \(\lambda\) optimo, tenemos lo siguiente:

bc <- boxcox(y ~ x)

Por lo tanto nuestro nuevo modelo:

library(ggplot2)
model_bc <- lm(y_transformed ~ x)
df <- data.frame(x = x, y_transformed = y_transformed)
ggplot(df, aes(x, y_transformed)) +
  geom_point() +
  geom_line(aes(y = predict(model_bc)), color = "red") +
  ggtitle("Ajuste lineal con transformación Box-Cox")

summary(lm_transformed)

## 
## Call:
## lm(formula = y_transformed ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.04491 -0.38752 -0.06783  0.37187  1.99376 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.96133    0.14298   6.724 1.19e-09 ***
## x            0.43473    0.02346  18.529  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6157 on 98 degrees of freedom
## Multiple R-squared:  0.7779, Adjusted R-squared:  0.7757 
## F-statistic: 343.3 on 1 and 98 DF,  p-value: < 2.2e-16

library(broom)
tidy(lm_transformed) %>% 
  mutate(
    p.value = format.pval(p.value, digits = 3),
    across(c(estimate, std.error, statistic), ~round(., 3))
  ) %>% 
  knitr::kable() %>% 
  kableExtra::kable_styling(bootstrap_options = "striped")

Table 1.3: Table 1.4: Summary del modelo transformado

term	estimate	std.error	statistic	p.value
(Intercept)	0.961	0.143	6.724	1.19e-09
x	0.435	0.023	18.529	< 2e-16

Así podemos notar que por (Tabla 1.3) ambos coeficientes son significantes para el modelo, además:

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
p1 <- ggplot(data.frame(x, y), aes(x, y)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  ggtitle("Datos Originales")

p2 <- ggplot(data.frame(x, y_transformed), aes(x, y_transformed)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, color = "blue") +
  ggtitle(paste("Datos Transformados (λ =", round(lambda, 2), ")"))

plot(lm_transformed, which = c(1))

Figure 1.5: Supuestos en la regresión datos transformados-1

plot(lm_transformed, which = c(2))

Figure 1.6: Supuestos en la regresión datos transformados-2

plot(lm_transformed, which = c(3))

Figure 1.7: Supuestos en la regresión datos transformados-3

ad_test <- ad.test(lm_transformed$residuals)
ad_test

## 
##  Anderson-Darling normality test
## 
## data:  lm_transformed$residuals
## A = 0.83135, p-value = 0.03093

Table 1.5: Resultados del test Anderson-Darling para los datos transformados

	Estadístico_AD	Valor_p
A	0.8313465	0.0309288

(Tabla 1.2)

No obstante vemos que por las Figuras (1.5,1.6,1.7) y en especial a la (Figura, 1.6), el supuesto de normalidad no es razonable, lo cual se confirma gracias a (Tabla 1.5)

Finalmente para el módelo multiple:

library(MASS)  
library(ggplot2)

set.seed(123)
x <- seq(1, 10, length.out = 100)
y <- exp(0.3 * x) + rgamma(100, shape = 2, rate = 0.5)  

# Ajuste cuadratico 
model_init <- lm(y ~ poly(x, degree = 2))  

# Transformación Box-Cox
bc <- boxcox(y ~ poly(x, degree = 2))  

lambda <- bc$x[which.max(bc$y)]
y_transformed <- (y^lambda - 1)/lambda
# Modelo transformado
model_bc <- lm(y_transformed ~ poly(x, degree = 2))  

# Gráfico del ajuste
ggplot(data.frame(x = x, y_transformed = y_transformed), aes(x, y_transformed)) +
  geom_point() +
  geom_line(aes(y = predict(model_bc)), color = "red", linewidth = 1) +
  ggtitle("Ajuste polinómico con transformación Box-Cox ") +
  theme_minimal()

Notemos lo siguiente:

shapiro.test(resid(model_bc))  # Test de normalidad

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model_bc)
## W = 0.95053, p-value = 0.0008995

qqnorm(resid(model_bc)); qqline(resid(model_bc))

2 Regresión no paramétrica

La regresión lineal ha servido como base fundamental para el desarrollo de métodos no paramétricos, en el enfoque lineal tradicional:

1. El usuario parte de un modelo tentativo para los datos
2. Verifica la razonabilidad de los supuestos:
   • Distribución normal de errores: \(\varepsilon_i \sim N(0, \sigma^2)\)
   • Linealidad de la relación
   • Homocedasticidad (varianza constante)
   • Independencia de los errores

En contraste, la regresión no paramétrica:

• Construye el modelo directamente desde los datos
• No requiere supuestos rígidos sobre la forma funcional

2.1 Métodos

El análisis de regresión no paramétrica comparte los objetivos de su contraparte paramétrica: estimar y contrastar propiedades de la función de regresión que describe la relación entre variables. Sin embargo:

\[ y = f(x) + \varepsilon \]

• La única condición sobre \(\varepsilon\) es que sean v.a.i.i.d con medía 0 y varianza \(\sigma^2\)

• La función de regresión \(f\) es desconocida y no se impone una forma funcional específica

• Se requieren supuestos mínimos sobre \(f\):

  • Continuidad (como condición básica)
  • Existencia de \(\lambda\) derivadas continuas (dependiendo del método)

Estas condiciones de suavidad permiten controlar el comportamiento local de las estimaciones.

2.2 Suavización

Supongamos que disponemos de una variable explicativa \(X\) y de una respuesta \(Y\) y que el modelo es aplicable al par \((X, Y )\), así sea f una función perteneciente a algún espacio de funciones \(\Omega\), con dicho espacio dependiendo de que tan suave es razonable asumir a \(f\),para estimar puntualmente a \(f\) en \(x\): \[\hat{f}(x)=\sum_{i=1}^nK(x,x_i;\lambda) y_i \] Es decir una suma ponderada donde \(K(x,x_i;\lambda)\) son pesos que dependen de \(\{x_i, i=1,..,n\}\) y del parámetro \(\lambda\) que es definido por el diseñador.

La siguiente imagen representa la probabilidad \(\mathbb{P}(x - h < X < x + h)\) y su aproximación mediante una densidad:

[{r mi-imagen, fig.cap="Aproximación de la densidad", echo=FALSE, out.width="70%", fig.align='center'} knitr::include_graphics("C:/Users/cyr/Documents/Aproxx.png")] : #

El estimador de densidad de núcleo se basa en una aproximación local vea ?? de la densidad \(f(x)\) utilizando los datos muestrales \(X_1, \dots, X_n \sim F\) (i.i.d.).

La idea fundamental es que la probabilidad de que una observación \(X\) caiga en un intervalo pequeño alrededor de un punto \(x\) puede aproximarse, para algún \(h > 0\) conocido como parámetro de suavizamiento o ancho de banda.

Para estimar esta probabilidad usando los datos disponibles, se puede contar la proporción de observaciones que caen dentro del intervalo \([x - h, x + h]\): Cuando el ancho de banda \(h\) es pequeño, podemos aproximar la integral usando el valor de \(f\) en el punto central \(x\):

\[ \int_{x-h}^{x+h} f(t) \, dt \approx 2h \cdot f(x) \] Pues por el Teorema del Valor Medio para integrales, existe \(c \in [x-h, x+h]\) tal que: \[ \int_{x-h}^{x+h} f(t) \, dt = f(c) \cdot 2h \] Cuando \(h \to 0\), \(c \to x\).

Usando la función de distribución empírica \(F_n\), estimamos la probabilidad como la proporción de puntos muestrales en el intervalo:

\[ \mathbb{P}(x - h \leq X \leq x + h) \approx \frac{1}{n} \sum_{i=1}^n \mathbf{1}_{\{|X_i - x| \leq h\}} \]

Igualando las expresiones:

\[ 2h \cdot f(x) \approx \frac{1}{n} \sum_{i=1}^n \mathbf{1}_{\{|X_i - x| \leq h\}} \]

Despejando \(f(x)\) obtenemos el estimador inicial:

\[ f(x) \approx \frac{1}{2h} \cdot \frac{1}{n} \sum_{i=1}^n \mathbf{1}_{\{|X_i - x| \leq h\}} = \hat{f}_h(x) \] Este estimador puede verse como una versión suavizada de un histograma centrado en \(x\). Para mejorar la suavidad, se reemplaza la función indicadora \(\mathbf{1}_{\{|X_i - x| \leq h\}}\) por alguna otra una función , típicamente una función simétrica y suave como la función gaussiana.

Así se obtiene el :

\[ \hat{f}_h(x) = \frac{1}{nh} \sum_{i=1}^n K\left( \frac{x - X_i}{h} \right). \]

Este estimador distribuye la “masa” de cada observación \(X_i\) a lo largo del eje usando el núcleo \(K\), centrado en \(X_i\) y escalado por \(h\), el resultado es una función suavizada que estima la densidad \(f(x)\).

Para el parámetro \(h\)

• \(h\) pequeño: Alta varianza (sobreajuste)
• \(h\) grande: Alto sesgo (subajuste)

Si los pesos \(K(x, xi; λ)\), \(i = 1, . . . , n\) son tales que:

• Si \(K(x, x_i; λ)\geq0\quad i=1,2,..,n\)
• Si \(\sum_{i=1}^n K(x, x_i; λ)=1\) Entonces el suavizador se dice normalizado

De lo contrario para los valores **dentro de la vecindad* se incluyen constante \(c_i\) tal que:

-\(\sum_{i=1}^n c_iK(x, x_i; λ)=1\)

Así notemos lo siguiente, para estimar una densidad, necesitamos de lo siguiente:

• \(n\) datos.
• \(h\) un ancho de banda

¿Cual es la configuración optima?

• Para 5000 datos y 2 disintias \(h\)

set.seed(123)
datos_unif <- runif(5000, min = -3, max = -1)

# Función de densidad teórica
densidad_teorica_unif <- function(x) {
  dunif(x, min = -3, max = -1)
}

# Kernel rectangular
kernel_rectangular <- function(x, centro, h) {
  ifelse(abs(x - centro) <= h, 1/(2*h), 0)
}

# Estimador de densidad
estimador_densidad <- function(x, datos, h) {
  mean(kernel_rectangular(x, datos, h))
}

# Rejilla de valores
rejilla <- seq(-4, 0, length.out = 500)

h1 <- 0.1
estimada_unif_h1 <- sapply(rejilla, estimador_densidad, datos = datos_unif, h = h1)

h2 <- 0.001
estimada_unif_h2 <- sapply(rejilla, estimador_densidad, datos = datos_unif, h = h2)

par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))

# Gráfico 1: h = 0.2
plot(rejilla, densidad_teorica_unif(rejilla), type = "l", lwd = 3, col = "red",
     main = paste("Kernel Rectangular (h =", h1, ")"),
     xlab = "x", ylab = "Densidad", ylim = c(0, 0.6))
lines(rejilla, estimada_unif_h1, lwd = 2, col = "blue")
rug(datos_unif, col = "gray50")
legend("topright", 
       legend = c("Teórica U(-3,-1)", paste("Kernel (h =", h1, ")")),
       col = c("red", "blue"), lwd = c(3, 2), cex = 0.8)

# Gráfico 2: h = 0.05
plot(rejilla, densidad_teorica_unif(rejilla), type = "l", lwd = 3, col = "red",
     main = paste("Kernel Rectangular (h =", h2, ")"),
     xlab = "x", ylab = "Densidad", ylim = c(0, 0.6))
lines(rejilla, estimada_unif_h2, lwd = 2, col = "blue")
rug(datos_unif, col = "gray50")
legend("topright", 
       legend = c("Teórica U(-3,-1)", paste("Kernel (h =", h2, ")")),
       col = c("red", "blue"), lwd = c(3, 2), cex = 0.8)

• Para 50 datos y 2 disintias \(h\)

set.seed(123)
datos_unif <- runif(50, min = -3, max = -1)

# Función de densidad teórica
densidad_teorica_unif <- function(x) {
  dunif(x, min = -3, max = -1)
}

# Kernel rectangular
kernel_rectangular <- function(x, centro, h) {
  ifelse(abs(x - centro) <= h, 1/(2*h), 0)
}

# Estimador de densidad
estimador_densidad <- function(x, datos, h) {
  mean(kernel_rectangular(x, datos, h))
}

# Rejilla de valores
rejilla <- seq(-4, 0, length.out = 500)

h1 <- 0.9
estimada_unif_h1 <- sapply(rejilla, estimador_densidad, datos = datos_unif, h = h1)

h2 <- 0.05
estimada_unif_h2 <- sapply(rejilla, estimador_densidad, datos = datos_unif, h = h2)

par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))

# Gráfico 1: h = 0.2
plot(rejilla, densidad_teorica_unif(rejilla), type = "l", lwd = 3, col = "red",
     main = paste("Kernel Rectangular (h =", h1, ")"),
     xlab = "x", ylab = "Densidad", ylim = c(0, 0.6))
lines(rejilla, estimada_unif_h1, lwd = 2, col = "blue")
rug(datos_unif, col = "gray50")
legend("topright", 
       legend = c("Teórica U(-3,-1)", paste("Kernel (h =", h1, ")")),
       col = c("red", "blue"), lwd = c(3, 2), cex = 0.8)

# Gráfico 2: h = 0.05
plot(rejilla, densidad_teorica_unif(rejilla), type = "l", lwd = 3, col = "red",
     main = paste("Kernel Rectangular (h =", h2, ")"),
     xlab = "x", ylab = "Densidad", ylim = c(0, 0.6))
lines(rejilla, estimada_unif_h2, lwd = 2, col = "blue")
rug(datos_unif, col = "gray50")
legend("topright", 
       legend = c("Teórica U(-3,-1)", paste("Kernel (h =", h2, ")")),
       col = c("red", "blue"), lwd = c(3, 2), cex = 0.8)

El riesgo cuadrático medio (\(R(\hat{f})\)) de un estimador \(\hat{f}\) se define como: \[ R(\hat{f}) = \mathbb{E}\left[(f(x) - \hat{f}(x))^2\right]. \]

Bajo condiciones de regularidad, puede descomponerse en: \[ R(\hat{f}) = \underbrace{\left(\mathbb{E}[\hat{f}(x)] - f(x)\right)^2}_{\text{Sesgo al cuadrado } (B^2)} + \underbrace{\text{Var}(\hat{f}(x))}_{\text{Varianza}}. \]

Demostración:

\[\begin{align*} R(\hat{f}) &= \mathbb{E}[(f(x) - \hat{f}(x))^2] \\ &= \mathbb{E}[(f(x) - \mathbb{E}[\hat{f}(x)] + \mathbb{E}[\hat{f}(x)] - \hat{f}(x))^2] \\ &= \underbrace{(f(x) - \mathbb{E}[\hat{f}(x)])^2}_{B^2} + \underbrace{\mathbb{E}[(\hat{f}(x) - \mathbb{E}[\hat{f}(x)])^2]}_{\text{Var}} + 2 {\mathbb{E}[(f(x) - \mathbb{E}[\hat{f}(x)])(\mathbb{E}[\hat{f}(x)] - \hat{f}(x))]}.\\ &= \underbrace{(f(x) - \mathbb{E}[\hat{f}(x)])^2}_{B^2} + \underbrace{\mathbb{E}[(\hat{f}(x) - \mathbb{E}[\hat{f}(x)])^2]}_{\text{Var}} \end{align*}\]

set.seed(123)
datos <- runif(50, -3, -1)
rejilla <- seq(-3.5, -0.5, length.out = 500)


kernel_rect <- function(x, c, h) ifelse(abs(x-c) <= h, 1/(2*h), 0)
estimador <- function(x, datos, h) rowMeans(sapply(datos, kernel_rect, x=x, h=h))

# riesgo 
riesgo <- function(h, datos, grid) {
  est <- estimador(grid, datos, h)
  teor <- dunif(grid, -3, -1)
  mse <- mean((est - teor)^2) + mean(apply(sapply(datos, kernel_rect, x=grid, h=h), 1, var))/length(datos)
  return(mse)
}

# Optimización
h_opt <- optimize(riesgo, c(0.05, 1.5), datos=datos, grid=rejilla)$minimum
cat("h óptimo:", round(h_opt, 3), "\n")

## h óptimo: 0.426

plot(h <- seq(0.05, 1.5, 0.02), sapply(h, riesgo, datos=datos, grid=rejilla), type='l',
     main="Riesgo vs h", xlab="h", ylab="Riesgo")
abline(v=h_opt, col=2, lty=2)

plot(rejilla, dunif(rejilla, -3, -1), type='l', col=2, lwd=2, ylim=c(0,0.6),
     main=paste("Estimación (h =", round(h_opt,3),")"))
lines(rejilla, estimador(rejilla, datos, h_opt), col=4)
rug(datos)

Confirmando el resultado, tenemos la Regla de Silverman \[ h_{opt}\approx n^{-1/5} \]

Que en nuestra caso sería, para \(50\) datos \[ h_{opt} \approx 0.4573\]

2.2.1 Propiedades kernel

Si ahora le pedimos a \(K: \mathbb{R} \to \mathbb{R}^+\) sea una función continua que cumple:

• Normalización: \(\int_{-\infty}^{\infty} K(u) \, du = 1\)
• Simetría: \(K(-u) = K(u)\)
• Colas : \(\lim_{|u|\to\infty} K(u) = 0\)

De esta manera se \(f_n(x)\) serán una densidad, pues: \[\begin{align*} \int_{-\infty}^{\infty}\hat{f_h} (x)&= \int_{-\infty}^{\infty} \frac{1}{nh} \sum_{i=1}^n K\left( \frac{x - X_i}{h} \right)dx\\ &=\dfrac{1}{nh}\sum_{i=1}^n\int_{-\infty}^{\infty}K\left( \frac{x - X_i}{h} \right)dx\\ &=\dfrac{1}{nh}\sum_{i=1}^n h\int_{-\infty}^{\infty}K\left(u\right)du\\ &=\dfrac{1}{nh}\sum_{i=1}^n h=1 \\ \end{align*}\]

Así en base a Berrendero (2023) para evaluar la calidad del estimador de densidad \(\hat{f}(x)\) en un punto fijo \(x\), se consideran dos medidas fundamentales:

1. Sesgo:
   \[ \text{Sesgo}[\hat{f}(x)] = E[\hat{f}(x)] - f(x) \]

2. Varianza:
   \[ \text{Var}[\hat{f}(x)] = E\left[(\hat{f}(x) - E[\hat{f}(x)])^2\right] \]

La consistencia débil del estimador \(\hat{f}(x)\) requiere que:

• El sesgo converja a cero: \(\text{Sesgo}[\hat{f}(x)] \to 0\)
• La varianza converja a cero: \(\text{Var}[\hat{f}(x)] \to 0\)

De esta manera tenemos que :

• Para reducir el sesgo: \(h_n \to 0\) cuando \(n \to \infty\)
• Para controlar la varianza: \(nh_n \to \infty\)

Así notemos lo siguiente ,el estimador kernel de densidad es:

\[ \hat{f}_h(x) = \frac{1}{nh}\sum_{i=1}^n K\left(\frac{x-X_i}{h}\right) \]

Como \(X_i\) son \(v.a.i.i.d\): \[\begin{align*} \mathbb{E}[\hat{f}_h(x)] &= \frac{1}{h}\mathbb{E}\left[K\left(\frac{x-X}{h}\right)\right] \\ &= \frac{1}{h}\int K\left(\frac{x-t}{h}\right)f(t)dt \\ &\text{Con el cambio de variable } u = \frac{x-t}{h}: \\ &= \int K(u)f(x-hu)du \end{align*}\]

Usando expansión de Taylor para \(h \to 0\):

\[\begin{align*} f(x-hu) &\approx f(x) - huf'(x) + \mathcal{O}(h^2) \\ \mathbb{E}[\hat{f}_h(x)] &\approx f(x)\underbrace{\int K(u)du}_{=1} - hf'(x)\underbrace{\int uK(u)du}_{=0} + \mathcal{O}(h^2) \\ &\to f(x) \text{ cuando } h \to 0 \end{align*}\]

Por otro lado la varianza, dado que \(X_i\) son \(v.a.i.i.d\)

\[\begin{align*} \text{Var}[\hat{f}_h(x)] &= \text{Var}\left[\frac{1}{nh}\sum_{i=1}^n K\left(\frac{x-X_i}{h}\right)\right] \\ &= \frac{1}{n^2h^2}\sum_{i=1}^n \text{Var}\left[K\left(\frac{x-X_i}{h}\right)\right] \quad \\ &= \frac{1}{nh^2}\text{Var}\left[K\left(\frac{x-X}{h}\right)\right] \end{align*}\]

Desarrollamos ahora la varianza del término kernel:

\[\begin{align*} \text{Var}\left[K\left(\frac{x-X}{h}\right)\right] &= \mathbb{E}\left[K^2\left(\frac{x-X}{h}\right)\right]_A - \left(\mathbb{E}\left[K\left(\frac{x-X}{h}\right)\right]\right)^2_B \\ &= \underbrace{\int K^2\left(\frac{x-t}{h}\right)f(t)dt}_{\text{A}} - \underbrace{\left(\int K\left(\frac{x-t}{h}\right)f(t)dt\right)^2}_{\text{ B}} \end{align*}\]

Desarrollo del Término A (cambio de variable u = (x-t)/h):

\[\begin{align*} \int K^2\left(\frac{x-t}{h}\right)f(t)dt &= h\int K^2(u)f(x-hu)du \\ &\approx hf(x)\int K^2(u)du \quad \text{(para h pequeño, f(x-hu) ≈ f(x))} \end{align*}\]

Desarrollo del Término B :

\[\begin{align*} \left(\int K\left(\frac{x-t}{h}\right)f(t)dt\right)^2 &= \left(h\int K(u)f(x-hu)du\right)^2 \\ &\approx \left(hf(x)\int K(u)du\right)^2 \\ &= h^2f^2(x) \end{align*}\] \[\begin{align*} \text{Var}[\hat{f}_h(x)] &\approx \frac{1}{nh^2}\left[hf(x)\int K^2(u)du - h^2f^2(x)\right] \\ &= \frac{f(x)}{nh}\int K^2(u)du - \frac{f^2(x)}{n} \end{align*}\]

Para h → 0 y nh → ∞, el primer término domina (es de orden O(1/(nh)) mientras que el segundo es O(1/n), por lo que:

\[\begin{align} \text{Var}[\hat{f}_h(x)] &\approx \frac{f(x)}{nh}\int K^2(u)du \quad \text{(aproximación asintótica)} \end{align}\]

1. \(h \to 0\) (reduce sesgo)
2. \(nh \to \infty\) (reduce varianza)

2.3 Estimador con series de Fourir

Para establecer un marco general, asumiremos que disponemos de observaciones de la variable de respuesta \(Y\) para \(n\) valores predeterminados de una variable independiente \(X\). Las \(n\) observaciones bivariadas disponibles, denotadas \((x_1, y_1), \ldots, (x_n, y_n)\), siguen el modelo

\[\begin{equation} \label{eq:regresion} y_i = f(x_i) + \epsilon_i, \quad i = 1, \ldots, n \tag{2.1} \end{equation}\]

donde \(\epsilon = (\epsilon_1, \ldots, \epsilon_n)^T\) es un vector de errores aleatorios no correlacionados que tienen media cero y varianza común \(\sigma^2\), \(f\) es una función de regresión desconocida y se satisface que \(0 \leq x_1 < \cdots < x_n \leq 1\).

Nuestro propósito es estimar \(f\), para lo cual buscaremos construir estimadores lineales que puedan escribirse en la siguiente forma general, que para un \(\lambda\) dado es una combinación lineal de las observaciones \(y_i\), donde \(K(\cdot, x_i; \lambda), \, i = 1, \ldots, n\) es una colección de funciones ponderadoras que dependen de los \(x_i\) y de un parámetro de suavización denotado\(\lambda\):

\[\begin{equation} f_\lambda(x)=\sum_{i=1}^nK(x, x_i; \lambda)y_i \end{equation}\]

Al final del siglo XVIII, Fourier propuso un método para aproximar funciones periódicas usando combinaciones lineales de funciones trigonométricas sencillas. Sin embargo, aunque el método propuesto por Fourier se refiere a funciones continuas, sus ideas se pueden aplicar mejor en el contexto más amplio de las funciones cuadrado integrables.

Definición: EL espacio \(L_2[0,1]\) se define como el espacio funcional de todas las funciones cuadrado integrables en [0, 1], es decir \[L_2[0,1]=\{f:[0,1] \rightarrow \mathbb{R}|\int_{0}^1 |f(x)|^2dx<\infty\}\]

Asumamos entonces que \(f \in L_2[0, 1]\) y que puede representarse usando la expansión: \[\begin{equation} f(x) = \sum_{j=1}^{\infty} \beta_j f_j(x) \tag{2.2} \end{equation}\]

donde los \(\beta_j\), \(j = 1, 2, \ldots\) son los , definidos en el Anexo 6.0.1, y las funciones \(f_j\), \(j = 1, 2, \ldots\) conforman una colección de funciones en \(L_2[0, 1]\). Diremos que la colección de funciones \(\{f_j\}_{j=1}^\infty\) utilizada para construir la expansión (2.2) es una (CONS), la definición de CONS se encuentra en el Anexo 6.0.1.

En consecuencia, el modelo (2.1) puede representarse como:

\[\begin{equation} y_i = \sum_{j=1}^{\infty} \beta_j f_j(x_i) + \varepsilon_i, \quad i = 1, 2, \ldots, n \end{equation}\]

Esto significa que los datos siguen un modelo lineal con una cantidad infinita de coeficientes de regresión desconocidos. Ahora bien, un resultado conocido sobre los coeficientes de Fourier, llamado la Relación de Parseval, definida como \[ \sum_{j=1}^{\infty} \beta_j^2 = \|f\|^2. \] Esta relación nos permite concluir que los coeficientes \(\beta_j\) deben decaer a cero eventualmente. Entonces uno podría entonces asumir que existe un entero \(\lambda\) tal que:

\[ f \approx \sum_{j=1}^{\lambda} \beta_j f_j \] y, por tanto, podríamos escribir la siguiente aproximación:

\[\begin{equation} y_i \approx \sum_{j=1}^{\lambda} \beta_j f_j(x_i) + \varepsilon_i, \quad i = 1, 2, \tag{2.3} \end{equation}\]

Este modelo (2.3) tiene la forma de un modelo lineal, por lo que una posible solución al problema de la estimación de \(f\) sería estimar los coeficientes de Fourier \(\{\beta_j\}_{j=1}^{\lambda}\) utilizando la técnica de mínimos cuadrados.

Para ello, definimos la matriz:

\[ X_\lambda = \{f_j(x_i)\}_{i=1,\ldots,n;\; j=1,\ldots,\lambda} \]

Si la matriz \((X_\lambda^\top X_\lambda)^{-1}\) existe, entonces el estimador de \(\boldsymbol{\beta}\) tendrá la forma general:

\[\begin{equation} \boldsymbol{\beta}_\lambda = (X_\lambda^\top X_\lambda)^{-1} X_\lambda^\top \mathbf{y} \tag{2.4} \end{equation}\]

con \(\mathbf{y} = (y_1, y_2, \ldots, y_n)^\top\). Entonces, nuestro estimador de \(f\) será:

\[\begin{equation} f_\lambda(x) = \sum_{j=1}^{\lambda} \beta_{\lambda j} f_j(x) \tag{2.5} \end{equation}\]

Supongamos que decidimos utilizar este estimador de \(f\). Quedarían en ese caso algunas dudas por resolver, recordando que \(\sum_{j=1}^{\lambda} \beta_j f_j\) es solamente una aproximación a \(f\) y que la verdadera función de regresión es, en realidad, de dimensión infinita.

2.3.0.1 Estimacion de la varianza

Una vez estimada la función de regresión, el análisis de regresión requiere estimar la varianza \(\sigma^2\) para poder continuar con el proceso de inferencia estadística. Un primer acercamiento podría ser adaptar el estimador de varianza del análisis de regresión lineal:

\[\begin{equation} \sigma^2_\lambda = \frac{\text{SCE}_\lambda}{n - \lambda} \tag{2.6} \end{equation}\]

donde la suma de cuadrados de los residuales \(\text{SCE}_\lambda\) se define como:

\[\begin{equation} \text{SCE}_\lambda = \sum_{i=1}^n \left(y_i - f_\lambda(x_i)\right)^2 \end{equation}\]

Sin embargo, esta definición tiene algunas dificultades, porque el estimador \(\sigma^2_\lambda\) definido en (2.6) depende del valor de \(\lambda\), que ya se ha utilizado para estimar \(f\). A diferencia del estimador de varianza en el modelo de regresión lineal, el estimador \(\sigma^2_\lambda\) no es, en general, insesgado. El sesgo del estimador \(\sigma^2_\lambda\) es:

\[\begin{equation} B(\sigma^2_\lambda) = \frac{f^\top (I - S_\lambda) f}{n - \lambda} \end{equation}\]

donde \(f = \big(f(x_1), f(x_2), \ldots, f(x_n)\big)^\top\) y la matriz de suavizamiento es:

\[ S_\lambda = X_\lambda (X_\lambda^\top X_\lambda)^{-1} X_\lambda^\top \]

Estudios previos concluyen que una elección de \(\lambda\) que minimice el sesgo del estimador \(\sigma^2_\lambda\) no necesariamente conduce a una buena selección del estimador \(f_\lambda\) de \(f\).

Una estimación de la varianza \(\sigma^2\) que no dependa de \(\lambda\) permitiría hacerla independiente del estimador de \(f\). En este contexto de los modelos no paramétricos se han propuesto varios estimadores de la varianza que no dependen de la elección de \(\lambda\). Algunos estimadores son:

• Estimador de Rice, fue propuesto por John Rice en 1984, y se define: \[\begin{equation} \sigma^2_R = \frac{1}{2(n - 1)} \sum_{i=2}^n (y_i - y_{i-1})^2 \end{equation}\] EL estimador \(\sigma^2_R\) depende en esencia de las contribuciones de los errores aleatorios que contienen toda la información acerca de \(\sigma^2\).
• Estimador de GSJS, Gasser, Sroka y Jennen-Steinmetz, apartir de la idea de Racer sugieren utilizar residuales obtenidos a partir de rectas ajustadas con tres puntos consecutivos para obtener otro estimador de la varianza. \[\begin{equation} \sigma^2_{\text{GSJS}} = \frac{1}{n - 2} \sum_{i=2}^{n-1} \tilde{\varepsilon}_i^2 \tag{2.7} \end{equation}\] Donde los \(\tilde{\varepsilon}_i\) se definen como los pseudo-residuales \[\begin{equation} \tilde{\varepsilon}_i = \frac{y_i - A_i y_{i-1} - B_i y_{i+1}}{\sqrt{1 + A_i^2 + B_i^2}}, \quad i = 2, \dots, n-1 \label{eq:residuales} \end{equation}\] con \[\begin{equation} A_i = \frac{x_{i+1} - x_i}{x_{i+1} - x_{i-1}}, \quad B_i = \frac{x_i - x_{i-1}}{x_{i+1} - x_{i-1}} \end{equation}\] El estimador GSJS busca resolver una limitación en los que se produzcan cambios muy fuertes de la función de regresión en subintervalos muy pequeños de \(X\). La solución en este caso consiste en medir la diferencia entre cada \(y_i\) y la recta que une sus dos vecinos más cercanos \(y_{i-1}\) y \(y_{i+1}\) y usar estos pseudo-residuales para obtener un estimador insesgado de \(\sigma^2\).
• Estimasor HKT Hall, Kay y Titterington proponen otro estimador basado en diferencias sucesivas. Sean \(m_1\) y \(m_2\) dos enteros no negativos tales que \(m_1 + m_2 = r\). El estimador \(\sigma^2_{\text{HKT}}\) se basa en sucesiones de diferencias de orden \(r\) definidas en el Anexo 6.0.1.3. El estimador \(\sigma^2_{\text{HKT}}\) se define en general como: \[\begin{equation} \sigma^2_{\text{HKT}} = \frac{1}{n - r} \sum_{l = m_1 + 1}^{n - m_2} \left(\Delta_l\right)^2=\frac{1}{n - r} \sum_{l = m_1 + 1}^{n - m_2} \left(\sum_{k = -m_1}^{m_2} d_k y_{k + l}\right)^2 \label{eq:hkt_general} \end{equation}\]

El estimador HKT es un estimador insesgado de \(\sigma^2\) que generaliza los dos estimadores previos de Rice y GSJS. El estimador de Rice coincide con el estimador HKT cuando \(r = 1\), con \(m_1 = 0\) y \(m_2 = 1\). Y el estimador GSJS se obtiene a partir del estimador HKT con \(m_1 = 0\) y \(m_2 = 2\), si el diseño es equidistante.

2.3.0.2 Algunas ideas sobre cómo conducir inferencias

Tal como en los modelos lineales, no parece razonable conservar una función \(f_j\) que tenga asociado un \(\beta_j = 0\). Así que una primera prueba formal del análisis de regresión será contrastar las hipótesis \[ H_0 : \beta_j = 0\quad\text{vs}\quad H_1 : \beta_j \neq 0 \] Un posible estadístico de prueba, adaptado de los modelos lineales, sería: \[ Z_j = \frac{\beta_j}{\sqrt{(X^\top_\lambda X_\lambda)_j \sigma^2}} \] donde \((X^\top_\lambda X_\lambda)_j\) es el elemento \(j\) de la diagonal principal de la matriz \(X^\top_\lambda X_\lambda\). En lugar de \(\sigma^2\) podríamos utilizar uno de los estimadores propuestos en la Sección 2.3.0.1 De igual forma, un intervalo de confianza para \(f\) sería: \[ f_\lambda(x) \pm Z_{\alpha/2} \sqrt{\sigma^2 f_\lambda(x)^\top (X^\top_\lambda X_\lambda)^{-1} f_\lambda(x)} \] con \(f_\lambda(x) = (f_1(x), \dots, f_\lambda(x))^\top\). El análisis de homocedasticidad y demás pruebas de diagnóstico se podrían tomar prestados de los métodos clásicos de análisis de regresión.

2.3.0.3 Ejemplo

El siguiente ejemplo fue tomado de Ochoa (2020b), los datos representados en la siguinete tabla corresponden a mediciones del promedio horario de Monóxido de Carbono (CO) en una estación ubicada en el centro de la ciudad de Cali, Colombia (Calle 15), el 1 de marzo de 2004.

Table 2.1: Mediciones de CO por hora

HORA	CO	HORA	CO
0	2.24	12	3.22
1	1.53	13	2.94
2	1.80	14	2.13
3	0.87	15	2.88
4	1.21	16	1.95
5	1.02	17	1.51
6	1.39	18	1.54
7	1.90	19	1.60
8	2.49	20	1.22
9	3.33	21	1.37
10	2.82	22	1.96
11	3.34	23	1.66

Nuestro propósito es ajustar un modelo del tipo (2.1) a este conjunto de observaciones de un día, para ello hagamos un reescalamiento de las horas para estar en el intervalo [0,1].

Supongamos que disponemos de la siguiente sucesión ortonormal completa (CONS) de funciones: \[\begin{align*} f_1(x) &= 1, \\ f_j(x) &= \sqrt{2} \cos((j - 1)\pi x), \quad j = 2, 3, \dots \end{align*}\]

Usando esto, nuestro estimador (2.5) tomaría la siguiente forma: \[\begin{equation} f_\lambda(x) = \sum_{j=1}^{\lambda} \beta_j^\lambda \sqrt{2} \cos((j - 1)\pi x) \end{equation}\]

Podemos estimar los \(\beta_{\lambda j}\) utilizando el método de mínimos cuadrados tal como se describió en (2.4), supongamos que para los datos decidimos estimar \(f\) usando \(\lambda=6\), es decir queremos aproximar \(f\) con las primeras 6 funciones de la CONS de cosenos:

\[\begin{equation} f_6(x) = \beta_{61} f_1 + \beta_{62} f_2 + \beta_{63} f_3 + \beta_{64} f_4+ \beta_{65} f_5+ \beta_{66} f_6 \end{equation}\]

# Datos
hora <- 0:23
co <- c(2.24, 1.53, 1.80, 0.87, 1.21, 1.02, 1.39, 1.90,
        2.49, 3.33, 2.82, 3.34, 3.22, 2.94, 2.13, 2.88,
        1.95, 1.51, 1.54, 1.60, 1.22, 1.37, 1.96, 1.66)

# Reescalar las horas al intervalo [0,1]
hora_reescalada <- (hora - min(hora)) / (max(hora) - min(hora))

# Definir las funciones f_j(x) para j = 2,...,6
f2 <- function(x) sqrt(2) * cos(1 * pi * x)
f3 <- function(x) sqrt(2) * cos(2 * pi * x)
f4 <- function(x) sqrt(2) * cos(3 * pi * x)
f5 <- function(x) sqrt(2) * cos(4 * pi * x)
f6 <- function(x) sqrt(2) * cos(5 * pi * x)

#  Matriz X_\lambda
matriz <- cbind(
  1,  # Columna de unos (intercepto)
  f2(hora_reescalada),
  f3(hora_reescalada),
  f4(hora_reescalada),
  f5(hora_reescalada),
  f6(hora_reescalada)
)

# Beta
beta <- solve(t(matriz) %*% matriz) %*% t(matriz) %*% co

# Mostrar los estimadores
print(beta)

##              [,1]
## [1,]  2.007396347
## [2,] -0.032968097
## [3,] -0.562119784
## [4,]  0.002013006
## [5,]  0.380031077
## [6,]  0.131421828

# Valores
x_nuevo <- seq(0, 1, length.out = 100)

#  Estimación de f
funcion_suma <- function(x) {
  beta[1] + 
    beta[2] * sqrt(2) * cos(1 * pi * x) + # f2
    beta[3] * sqrt(2) * cos(2 * pi * x) + # f3
    beta[4] * sqrt(2) * cos(3 * pi * x) + # f4
    beta[5] * sqrt(2) * cos(4 * pi * x) + # f5
    beta[6] * sqrt(2) * cos(5 * pi * x)   # f6
}

# Evaluar en los nuevos puntos
y_ajustado <- sapply(x_nuevo, funcion_suma)

# Graficar datos con la función ajustada
plot(hora_reescalada, co, 
     main = "Ajuste de funciones base coseno",
     xlab = "Hora (reescalada a [0,1])",
     ylab = "Concentración de CO",
     pch = 19, col = "blue", cex = 0.7)

lines(x_nuevo, y_ajustado, col = "red", lwd = 2)

2.3.1 El modelo de regresión polinómica

Asumamos por un momento que disponemos de una función de regresión \(f\) conocida. De acuerdo con el Teorema de Taylor, si tenemos una función continua f en el intervalo cerrado \([a, b]\) que tiene \(\lambda+1\) derivadas continuas en el intervalo (a, b), entonces si \(x\) y \(c\) son dos puntos en (a, b), se sigue que \(f(x)\) se puede representar como: \[\begin{equation} f(x)=f(c)+(x-c)f'(c)+\frac{(x-c)^2}{2!}f''(c)+\dots+ \frac{(x-c)^\lambda}{\lambda!}f''(c)+r_\lambda(\xi). \tag{2.8} \end{equation}\]

Donde el término \(r_\lambda(\xi)\) depende de la derivada \(f^{\lambda+1}\) de \(f\) para cierto \(\xi\) entre \(x\) y \(c\), y se conoce como el \(\lambda\)-residuo. La expresión de \(f\) usando el Teorema de Taylor @ref{eq:Taylor} puede reescribirse como:

\[\begin{equation} f(x)=\sum_{j=0}^\lambda\theta_jf_j(x)+r_\lambda(x). \tag{2.9} \end{equation}\]

donde \[ \theta_1 = f(c), \quad \theta_j = \frac{f^{(j-1)}(c)}{(j-1)!}, \quad j = 2, \ldots, \lambda, \quad \text{y}\quad f_j(x) = (x - c)^{j-1}, \quad j = 1, \ldots, \lambda. \]

Es decir, la función \(f\) se puede representar como la suma de un polinomio \(\pi_\lambda(x)\) de grado \(\lambda\) más un residuo \(r_\lambda(x)\).

Ahora bien, en la medida que \(r_\lambda(x)\) sea más pequeño, entonces \(\pi_\lambda(x)\) se parecerá más a \(f(x)\). En resumen, existe un \(\lambda\) tal que: \[ f(x) \approx \sum_{j=0}^{\lambda} \theta_j f_j(x) \]

Supongamos ahora que hemos colectado información sobre dos variables \(X\) y \(Y\) que se relacionan según el modelo (2.1), por lo que la función de regresión no es conocida. Si asumimos que \(f\) es una función de regresión que tiene \(\lambda + 1\) derivadas continuas y puede representarse en la forma (2.9), entonces, si \(r_\lambda(x)\) es uniformemente pequeño, podríamos escribir:

\[\begin{equation} f(x) \approx \sum_{j=0}^{\lambda} \theta_j f_j(x) + \varepsilon_j \tag{2.10} \end{equation}\]

Diremos entonces que nuestras observaciones se comportan aproximadamente según un modelo de regresión (no-paramétrica) polinómica, el cual presupone que los residuos \(r_\lambda(x_1), \ldots, r_\lambda(x_n)\) de la expansión de \(f\) en series de Taylor han sido incorporados en los errores aleatorios del modelo (2.10).

Para encontrar los coeficientes \(\theta_j\), se podría utilizar el método de mínimos cuadrados, o cualquier otro.

Por otra parte, si \(f\) se puede aproximar bien utilizando el polinomio \(\pi_\lambda(x)\), y si los residuos en los \(x_i\) son pequeños en comparación con los errores \(\varepsilon_i\), entonces el modelo (2.10) debería funcionar bien. En otro caso, podríamos estar en problemas con este modelo, porque como tanto los residuos como los errores son desconocidos.

3 Estimadores Kernel

3.1 Introducción

Para establecer un marco general, se asume que se dispone de observaciones de la variable de respuesta \(Y\) para \(n\) valores predeterminados de una variable independiente \(X\). Las \(n\) observaciones bivariadas disponibles, denotadas \((x_1, y_1), \dots, (x_n, y_n)\), siguen el modelo

\[ y_i = f(x_i) + \varepsilon_i, \quad i = 1, \dots, n \tag{3.1} \] donde \(\varepsilon = (\varepsilon_1, \dots, \varepsilon_n)^T\) es un vector de errores aleatorios no correlacionados que tienen media cero y varianza común \(\sigma^2\), \(f\) es una función de regresión desconocida y se satisface que \(0 \le x_1 < \cdots < x_n \le 1\).

Según (Ochoa 2020b), para efectos de presentación de los resultados, se asume que los valores de \(X\) se han elegidos así: \[ x_i = \frac{2i - 1}{2n}, \quad i = 1, 2, \dots, n \tag{3.2} \] El propósito es estimar \(f\) en (3.1), para lo cual se busca construir estimadores lineales que puedan escribirse en la siguiente forma general, que para un \(\lambda\) dado es una combinación lineal de las observaciones \(y_i\), donde \(K(\cdot; x_i; \lambda)\), \(i = 1, \ldots, n\) es una colección de funciones ponderadoras que dependen de los \(x_i\) y de un parámetro de suavización denotado \(\lambda:\)

\[ f_\lambda(x) = \sum_{i=1}^n K(x, x_i; \lambda) y_i \tag{3.3} \]

3.2 Estimadores Kernel

De acuerdo con Ochoa (2020b), se discute sobre los estimadores de series, donde se encontra que el estimador de cosenos puede representarse como un estimador lineal, que se construye con una colección de funciones de pesos, una para cada \(x_i\), que depende de \(\lambda:\)

\[ f_\lambda(x) = \frac{1}{n} \sum_{i=1}^{n} y_i K_\lambda(x; x_i) \tag{3.4} \] donde:

\[ K_\lambda(x; s) = \left[1 + \cos\left(\frac{\lambda \pi (x - s)}{2}\right)\right] D\left(\frac{x - s}{2}; \lambda - 1\right) + \left[1 + \cos\left(\frac{\lambda \pi (x + s)}{2}\right)\right] D\left(\frac{x + s}{2}; \lambda - 1\right) \tag{3.5} \] En (3.5) la función \(D\) se conoce como kernel de Dirichlet, que proviene del análisis de series de Fourier y se define como:

La función de asignación de pesos (3.5) que se utiliza en ese caso tiene la limitación de oscilar a medida que se aleja del punto de estimación \(x\), de tal manera que asigna pesos bajos a algunos puntos cercanos al punto de estimación y pesos más altos a puntos más alejados de este punto de estimación.

\[ D(u; k) = \frac{\sin(k \pi u)}{\sin(\pi u)} \]

Con este enfoque, la observación \(y_i\) asociada con un \(x_i\) cercano a \(x\) tendrá asignado un peso más alto, contribuyendo lo máximo posible a la estimación de la función de regresión en \(x\). Al mismo tiempo, observaciones \(y_i\) asociadas con \(x_i\) alejadas de \(x\) tendrán pesos cada vez más pequeños a medida que estén más alejadas de \(x\).

Una opción de construcción de estimadores que eviten esta limitación sería utilizar la misma idea pero con una función ponderadora que asigne consistentemente pesos más bajos a observaciones más alejadas del punto de estimación, al mismo tiempo que asigna pesos altos a las observaciones cercanas al punto deestimación. Una manera de producir estos estimadores es utilizar colecciones de funciones ponderadoras que tengan la forma general:

\[ K(x, x_i; \lambda) = \frac{1}{\lambda} K\left(\frac{x - x_i}{\lambda}\right), \quad i = 1, \dots, n \tag{3.6} \] donde \(K\) es una función simétrica con soporte en \([-1, 1]\) y su máximo en cero. Al parámetro \(\lambda\), que en los estimadores de series era el número de términos de la serie que se conservaban en la estimación de \(f\), lo llamaremos indistintamente parámetro de suavización o ancho de banda y cumplirá la misma función: servir de base para elegir el mejor estimador de \(f.\) Pero a diferencia de los estimadores de series, en los estimadores kernel \(\lambda\) no tiene que ser un entero, aunque sí se restringirá su valor a que sea cualquier real no negativo.

Con este modo de hacer, la observación \(y_i\) asociada con un \(x_i\) cercano a \(x\) tendrá asignado un peso más alto, contribuyendo de esta manera lo máximo posible a la estimación de la función de regresión en \(x\). Al mismo tiempo, observaciones \(y_i\) asociadas con \(x_i\) alejadas de \(x\) tendrán pesos cada vez más pequeños a medida que estén más alejados de \(x\). Para garantizar lo anterior, le exigiremos por lo pronto a la función \(K\) que tenga las siguientes dos propiedades:

\[ \int_{-1}^1 K(u) \, du = 1 \\ \int_{-1}^1 u K(u) \, du = 0 \tag{3.7} \]

La primera propiedad en (3.7) ayuda a garantizar que la suma de los pesos sea igual a 1; la segunda, a que \(K\) sea simétrica alrededor de cero.

A estas funciones se les llama funciones kernel. Y a los estimadores basados en estas funciones se les llama estimadores kernel.

Dos funciones kernel de uso común, que se conocen como:

• Kernel cuadrático (Epanechnikov): \[ K(u) = \frac{3}{4}(1 - u^2) \quad \text{si } |u| \le 1, \text{ y } 0 \text{ en otro caso} \]

• Kernel biponderado (biweight): \[ K(u) = \frac{15}{16}(1 - u^2)^2 \quad \text{si } |u| \le 1, \text{ y } 0 \text{ en otro caso} \]

# Definimos las funciones kernel
epanechnikov <- function(u) ifelse(abs(u) <= 1, 0.75 * (1 - u^2), 0)
biweight <- function(u) ifelse(abs(u) <= 1, (15/16) * (1 - u^2)^2, 0)
#gaussian <- function(u) dnorm(u)  # sin recorte

# Valores para graficar
u_vals <- seq(-1.5, 1.5, length.out = 500)

# Graficamos
plot(u_vals, epanechnikov(u_vals), type = "l", col = "blue", lwd = 2,
     ylim = c(0, 1), ylab = "K(u)", xlab = "u", main = "Funciones Kernel")
lines(u_vals, biweight(u_vals), col = "red", lwd = 2)
#lines(u_vals, gaussian(u_vals), col = "green", lwd = 2, lty = 2)
legend("topright", legend = c("Epanechnikov", "Biweight"),
       col = c("blue", "red"), lty = c(1, 1, 2), lwd = 2)

3.3 Estimador kernel de Priestley-Chao

Con el estimador de cosenos se concluyo que una forma particular del estimador es:

\[ f_{\lambda}(x) = \frac{1}{n} \sum_{i=1}^n K(x; x_i; \lambda)y_i \]

Si se incorpora ahora el kernel genérico, se obtiene el estimador kernel:

\[ f_{\lambda}(x) = \frac{1}{n \lambda} \sum_{i=1}^n K\left(\frac{x - x_i}{\lambda}\right)y_i \]

Este estimador kernel es un estimador kernel genérico que luce tal como el estimador de cosenos (3.5). Aunque ambos asignan pesos altos a puntos cercanos a \(x\), el kernel genérico reduce el peso más rápidamente con la distancia.

Este estimador fue propuesto por Priestley and Chao (1972) para el caso de diseños igualmente espaciados como el del diseño. Sin embargo, su desempeño no es muy bueno, como se ve en la Figura 3.1, donde se destaca la escasa habilidad del estimador en los extremos del diseño.

# Datos de hora y CO
hora <- 0:23
co <- c(2.24, 1.53, 1.80, 0.87, 1.21, 1.02, 1.39, 1.90, 2.49, 3.33, 2.82, 3.34,
        3.22, 2.94, 2.13, 2.88, 1.95, 1.51, 1.54, 1.60, 1.22, 1.37, 1.96, 1.66)

# Normalizar la hora al intervalo [0,1]
hora_cod <- hora / max(hora)

# Kernel cuadrático (biweight)
k_biweight <- function(u) {
  ifelse(abs(u) <= 1, (15/16) * (1 - u^2)^2, 0)
}

# Estimación kernel para un punto x dado
f_lambda <- function(x, x_i, y_i, lambda) {
  n <- length(x_i)
  sum(k_biweight((x - x_i)/lambda) * y_i) / (n * lambda)
}

# Evaluar la estimación en una malla de puntos
lambda <- 0.15
x_grid <- seq(0, 1, length.out = 200)
f_est <- sapply(x_grid, function(x) f_lambda(x, hora_cod, co, lambda))

# Graficar
plot(hora_cod, co,main="Estimación Kernel de Concentraciones de CO en Cali" ,pch = 19,col="blue" ,xlab = "HORA codificada", ylab = "CO")
lines(x_grid, f_est, type = "l")

Figure 3.1: Estimacion kernel del comportamiento diario del CO horario en la calle 15 de Cali usando el Kernel generico con el kernel cuadratico y con \(lambda= 0.15\)

El ancho de banda se elige de manera similar, basado en la minimización del criterio de validación cruzada generalizada GCV.

Si uno dispusiera de un segundo conjunto de \(n\) datos, estimaría \(f_\lambda\) con los primeros \(n\) datos y buscaría el \(\lambda\) que minimice \(P(\lambda)\) usando el segundo conjunto de datos. En la práctica, sin embargo, sería mejor estimar \(f\) utilizando los \(2n\) datos, lo que nos lleva a continuar nuestra búsqueda de un nuevo procedimiento para elegir \(\lambda\) que no dependa de \(\sigma^2\).

Una opción es dividir el conjunto de \(n\) observaciones en \(n\) sub-muestras de tamaño \(n - 1\) mediante el mecanismo de dejar por fuera una observación diferente cada vez. Si denotamos \(\hat{f}_\lambda^{(i)}\) a la estimación de \(f_i\) obtenida al suprimir de la muestra la observación \(i\), entonces la observación \(y_i\) sería una observación adicional que podríamos utilizar para construir un estimador de \(P(\lambda)\) que denotaremos \(\text{CV}(\lambda)\) y que llamaremos criterio de validación cruzada:

\[ \text{CV}(\lambda) = \frac{1}{n} \sum_{i=1}^n (y_i - \hat{f}_\lambda^{(i)})^2 \]

Calcular el criterio CV es complejo, porque requiere una gran carga de procesamiento computacional. Green and Silverman (1993) sugieren simplificar este cálculo utilizando en cambio:

\[ \text{CV}_{GS}(\lambda) = \frac{1}{n} \sum_{i=1}^n \left( \frac{y_i - \hat{f}_{\lambda, i}}{1 - S_{\lambda, ii}} \right)^2 \]

donde \(S_{\lambda, ii}\) es el elemento \(i\) de la matriz \(S_\lambda\).

Otra posible solución se debe a Wahba (1990), quien sugiere
utilizar otro criterio llamado validación cruzada generalizada
\(\text{GCV}(\lambda)\), definido por como:

\[ \text{GCV}_{GS}(\lambda) = \frac{n^{-1} \sum_{i=1}^n (y_i - \hat{f}_{\lambda, i})^2}{(1 - n^{-1} \text{tr}(S_\lambda))^2} \]

Comparando las expresiones se ve que los residuales obtenidos al eliminar una de las observaciones se pueden obtener a partir de los residuales mediante el mecanismo de dividirlos entre el factor \(1 - S_{\lambda, ii}\). La idea del criterio GCV es reemplazar estos factores por su promedio \(1 - n^{-1} \text{tr}(S_\lambda)\).

Como en el numerador aparece la suma de cuadrados de los residuales \(\text{RSS}\), una expresión más sencilla sería Eubank (1999):

\[ \text{GCV}(\lambda) = \frac{n^{-1} \text{RSS}(\lambda)}{(n^{-1} \text{tr}(I - S_\lambda))^2} \tag{3.8} \]

Aunque la palabra “generalizada” deja la impresión de que el segundo criterio generaliza el primero, esto no es en general cierto y se trata de criterios diferentes que permiten estimar el riesgo de predicción \(P(\lambda)\). Wahba (1990) justifica el uso de este criterio como un buen método de selección de \(\lambda\) demostrando que si \(n^{-1} \text{tr}(S_\lambda) < 1\), entonces la diferencia entre \(\mathbb{E}[\text{GCV}(\lambda)]\) y \(P(\lambda)\), relativa al tamaño de \(R(\lambda)\), será pequeña, especialmente para tamaños de muestra grande.

En este caso el estimador lucirá como \(\hat{f}_\lambda = S_\lambda y\) , donde un término \(S_{ij}\) de la matriz \(S_\lambda\) lucirá como:

\[ S_{ij} = \frac{1}{n \lambda} K\left( \frac{x_j - x_i}{\lambda} \right), \quad i, j = 1, \dots, n \]

Otro aspecto importante es la selección del kernel. Una alternativa sería identificar el kernel que minimiza el riesgo de \(\hat{f}_h\), una vez elegido el ancho de banda. Estudios previos (Benedetti 1975, Gasser & Müller 1979) identifican el kernel cuadrático de Epanechnikov (1969), adaptado para que sea positivo en \([-1, 1]\), como una buena elección.

Otros kernels propuestos incluyen:

• Kernel Gaussiano:

\[ K(u) = \frac{1}{\sqrt{2\pi}} e^{-u^2}, \quad u \in (-\infty, \infty) \]

• Versión simplificada del kernel Gaussiano:

\[ K(u) = e^{-u^2}, \quad u \in (-\infty, \infty) \]

• Kernel triponderado (triweight):

\[ K(u) = \frac{35}{32}(1 - u^2)^3, \quad u \in [-1, 1] \]

Sin embargo, la elección del kernel tiene baja influencia comparada con la elección del ancho de banda.

3.4 Estimador de Nadaraya-Watson

Varios intentos posteriores tratan de mejorar la habilidad del estimador genérico. El primero fue propuesto por Benedetti (1975), y también por Nadaraya (1964) y Watson (1964), para el caso en el que \(X\) es una variable aleatoria. Este estimador de conoce como estimador de Nadaraya-Watson y se define así:

\[ f_{\lambda}(x) = \frac{\sum_{i=1}^n K\left(\lambda^{-1}(x - x_i)\right)y_i}{\sum_{j=1}^n K\left(\lambda^{-1}(x - x_j)\right)} \]

Este estimador ajusta los pesos del kernel para que sumen 1, convirtiéndolo en un promedio ponderado.

# Normalizar la hora a [0, 1]
hora_cod <- hora / max(hora)

# Estimación kernel usando ksmooth (Nadaraya-Watson con kernel normal)
ks <- ksmooth(x = hora_cod, y = co, kernel = "normal", bandwidth = 0.15)

# Graficar resultados
plot(hora_cod, co, pch = 1, xlab = "HORA codificada", ylab = "CO")
lines(ks, col = "blue", lwd = 2)

La versión original parte de:

\[ f(x) = E(Y \mid X = x) \]

Requiriendo la distribución condicional \(h(y|x)\), y usando:

\[ f(x) = \int y h(y|x) \, dy = \int y \frac{h(x, y)}{h_X(x)} \, dy \tag{3.9} \]

Una solución es estimar densidades vía suavización. Härdle (1990b) y Simonoff (1996) discuten extensamente esta metodología.

Un estimador kernel-producto de \(h(x, y)\) es:

\[ \hat{h}(x, y) = \frac{1}{h_x h_y} \sum_{i=1}^n K_x\left(\frac{x - x_i}{h_x}\right) K_y\left(\frac{y - y_i}{h_y}\right) \]

Y un estimador de \(h_X(x)\) es:

\[\hat{h}(x) = \frac{1}{h_x} \sum_{i=1}^n K_x\left(\frac{x - x_i}{h_x}\right) \] Con estas ecuaciones se obtiene el estimador de Nadaraya-Watson.

3.5 Estimador de Gasser-Müller

Otro estimador kernel fue propuesto por Gasser & Müller (1979):

\[ \hat{f}(x) = \sum_{i=1}^n \left[ h^{-1} \int_{s_{i-1}}^{s_i} K\left(h^{-1}(x - s)\right) ds \right] y_i \]

con:

\[ s_0 = 0, \quad s_{i-1} \leq x_i \leq s_i, \quad i = 1, \dots, n - 1, \quad s_n = 1 \]

# Normalizar hora
x <- hora / max(hora)
y <- co
n <- length(x)
lambda <- 0.15

# Kernel normal
K <- function(u) dnorm(u)

# Crear los límites s_i
s <- numeric(n + 1)
s[1] <- x[1] - (x[2] - x[1]) / 2
s[n + 1] <- x[n] + (x[n] - x[n - 1]) / 2
for (i in 2:n) {
  s[i] <- (x[i - 1] + x[i]) / 2
}

# Estimador de Gasser-Müller
gasser_muller <- function(x_eval) {
  sum(sapply(1:n, function(i) {
    integral <- integrate(function(s_var) K((x_eval - s_var)/lambda), lower = s[i], upper = s[i + 1])$value
    (1 / lambda) * integral * y[i]
  }))
}

# Evaluar en una malla de puntos
x_grid <- seq(0, 1, length.out = 200)
f_gm <- sapply(x_grid, gasser_muller)

# Graficar
plot(x, y, pch = 1, xlab = "HORA codificada", ylab = "CO")
lines(x_grid, f_gm, type = "l", col = "red", lwd = 2)

Tenemos entonces tres estimadores kernel:

• Estimador de Priestley-Chao
• Estimador de Nadaraya-Watson
• Estimador de Gasser-Müller

Estos estimadores tienen propiedades asintóticas similares. asintóticas, por lo que en lo sucesivo no nos detendremos a identificar cuál de ellos está en uso. desigualmente espaciados, mientras que el estimador de Gasser-Muller está diseñado específicamente para este tipo de diseños desigualmente espaciados. El estimador de Nadaraya-Watson fue pensado por sus autores para diseños aleatorios y ha tomado la delantera entre la mayoría de los usuarios de los métodos de regresión no paramétrica por su flexibilidad y eficiencia.

3.6 Estimadores lineales localmente

Otro estimador común, que también es un estimador kernel, es el estimador de regresión local Cleveland (1979). Se sabe por Eubank (1999) que el estimador de Nadaraya-Watson coincide con la solución para \(\beta_0\) del problema de mínimos cuadrados ponderados:

\[ \sum_{i=1}^n (y_i - \beta_0)^2 K\left(\frac{x - x_i}{h}\right), \tag{3.9} \] lo que sugiere el ajuste de polinomios locales de orden \(p\), por ejemplo uno que minimice:

\[ \sum_{i=1}^n \left(y_i - \beta_0 - \dots - \beta_p(x - x_i)^p\right)^2 K\left(\frac{x - x_i}{h}\right). \tag{3.10} \] Al estimador resultante de minimizar la suma (3.10) se le conoce como estimador polinomial de regresión local de orden \(p\).

3.6.1 Implementación del estimador LOESS

Ilustraremos el caso de la solución del sistema (3.9), popularizado como estimador LOESS, que utiliza los \(k\) vecinos más cercanos:

1. Identificar los \(k\) vecinos más cercanos de \(x\) (denotados como \(N(x)\)).
2. Calcular \(\Delta(x) = \max_{x_i \in N(x)} |x - x_i|\) (distancia al vecino más alejado en \(N(x)\)).
3. Asignar pesos \(w_i\) usando el tri-cubo: \[ W\left(\frac{|x - x_i|}{\Delta(x)}\right), \quad \text{donde} \quad W(u) = \begin{cases} (1 - u^3)^3, & 0 \leq u \leq 1 \\ 0, & \text{otro caso} \end{cases} \]
4. Ajustar una recta por mínimos cuadrados ponderados de \(Y\) sobre \(X\) (usando los pesos \(w_i\)).

La estimación LOESS de \(f(x)\) toma el valor del intercepto (\(\hat{\beta}_0\)) de la recta ajustada \(\hat{\beta}_0 + \hat{\beta}_1 x\).

# Datos
hora <- 0:23
co <- c(2.24, 1.53, 1.80, 0.87, 1.21, 1.02, 1.39, 1.90, 2.49, 3.33, 2.82, 3.34,
        3.22, 2.94, 2.13, 2.88, 1.95, 1.51, 1.54, 1.60, 1.22, 1.37, 1.96, 1.66)

# Normalizar la hora a [0,1]
hora_cod <- hora / max(hora)

# Estimación LOESS (estimación lineal local con kernel Gaussiano)
modelo_loess <- loess(co ~ hora_cod, span = 0.5, degree = 2, family = "gaussian")

# Crear una malla para evaluar el modelo
x_grid <- seq(0, 1, length.out = 200)
f_loess <- predict(modelo_loess, newdata = data.frame(hora_cod = x_grid))

# Graficar
plot(hora_cod, co, pch = 1, xlab = "HORA codificada", ylab = "CO")
lines(x_grid, f_loess, col = "darkgreen", lwd = 2)

# Datos
hora <- 0:23
co <- c(2.24, 1.53, 1.80, 0.87, 1.21, 1.02, 1.39, 1.90, 2.49, 3.33, 2.82, 3.34,
        3.22, 2.94, 2.13, 2.88, 1.95, 1.51, 1.54, 1.60, 1.22, 1.37, 1.96, 1.66)

# Normalizar la hora a [0,1]
hora_cod <- hora / max(hora)

# Estimación LOESS (estimación lineal local con kernel Gaussiano)
modelo_loess <- loess(co ~ hora_cod, span = 0.2, degree = 2, family = "gaussian")

# Crear una malla para evaluar el modelo
x_grid <- seq(0, 1, length.out = 200)
f_loess <- predict(modelo_loess, newdata = data.frame(hora_cod = x_grid))

# Graficar
plot(hora_cod, co, pch = 1, xlab = "HORA codificada", ylab = "CO")
lines(x_grid, f_loess, col = "darkgreen", lwd = 2)

# Datos
hora <- 0:23
co <- c(2.24, 1.53, 1.80, 0.87, 1.21, 1.02, 1.39, 1.90, 2.49, 3.33, 2.82, 3.34,
        3.22, 2.94, 2.13, 2.88, 1.95, 1.51, 1.54, 1.60, 1.22, 1.37, 1.96, 1.66)

# Normalizar la hora a [0,1]
hora_cod <- hora / max(hora)

# Estimación LOESS (estimación lineal local con kernel Gaussiano)
modelo_loess <- loess(co ~ hora_cod, span = 0.5, degree = 1, family = "gaussian")

# Crear una malla para evaluar el modelo
x_grid <- seq(0, 1, length.out = 200)
f_loess <- predict(modelo_loess, newdata = data.frame(hora_cod = x_grid))

# Graficar
plot(hora_cod, co, pch = 1, xlab = "HORA codificada", ylab = "CO")
lines(x_grid, f_loess, col = "darkgreen", lwd = 2)

3.7 Estimación kernel multivariante

Todos los estimadores kernel discutidos en esta sección pueden extenderse a la estimación de funciones \(p\)-variadas. Supongamos que \(\mathbf{x} = (x_1, x_2, \dots, x_p)^T\). La solución es definir kernels multivariados como productos de kernels univariados:

\[ K(u_1, \dots, u_p) = \prod_{j=1}^p K_j(u_j), \tag{3.11} \]

El estimador de Priestley-Chao multivariado toma la forma:

\[ \hat{f}(\mathbf{x}) = \frac{1}{n} \sum_{i=1}^n y_i \prod_{h=1}^p h_h^{-1} K_h\left(h_h^{-1}(x_h - x_{ih})\right) \quad (3.16) \]

3.7.1 Desafíos en estimación multivariante

1. Complejidad computacional:
   • Requiere seleccionar múltiples parámetros de suavización (\(h_1, \dots, h_p\)).
   • Costo computacional crece exponencialmente con \(p\) (problema de curse of dimensionality).
2. Visualización:
   • Gráficos 3D son posibles para \(p=2\), pero no para \(p \geq 3\).
3. Maldición de la dimensionalidad:
   • En espacios de alta dimensión, los datos se vuelven escasos.

3.7.2 Forma general del estimador kernel multivariante

\[ \hat{f}_L(\mathbf{x}) = \frac{1}{n |L|} \sum_{i=1}^n K_p\left[L^{-1}(\mathbf{x} - \mathbf{x}_i)\right] \]

donde: - \(K_p: \mathbb{R}^p \to \mathbb{R}\) es un kernel multivariado (usualmente producto de kernels univariados). - \(L\) es una matriz \(p \times p\) de anchos de banda. - \(|L|\) es el determinante de \(L\). no-singular; y \(|L|\) es el valor absoluto del determinante de la matriz \(L\). La técnica más común para obtener \(K_p\) a partir de una función kernel univariada \(K\) es utilizar un producto de la forma (3.11).

Un ejemplo de estimación no paramétrica bivariada se presenta a continuación, en la que se estima el comportamiento del contaminante Ozono (O₃) dependiendo de la temperatura y la humedad relativa.

En este problema de la contaminación por ozono troposférico en Londres enero del 2002 , hemos ajustado un modelo bivariado de la forma:

\[ Y = f(X_1, X_2) + \varepsilon \]

donde \(Y\) representa el promedio diario de Ozono y \(X_1\) y \(X_2\) los promedios diarios de temperatura y humedad, respectivamente.

Figure 3.2: 3DPLOT

3.8 Inferencia en la estimación kernel

1. Primera condición
   Garantiza que el primer término de \(\mathbb{E}[\hat{f}_h(x)]\) en (3.22) es exactamente \(f(x)\), asegurando insesgabilidad asintótica en el interior del diseño.

2. Segunda condición
   Suprime el término de orden \(h\) en el sesgo (ver Ejercicio 2 al final del capítulo), reduciéndolo a \(\mathcal{O}(h^2 + n^{-1})\).
   Demostración:
   \[\int_{-1}^1 u K(u) du = 0 \implies M_1(x) = 0\]

3. Primera condición

   • \(M_2 \neq 0\) implica que \(K\) captura curvatura (\(f \in C^2[0,1]\)).
     
   • El sesgo decae óptimamente como \(\mathcal{O}(h^2)\) para funciones suaves.
     

4. Segunda condición
   \(V < \infty\) garantiza varianza finita:
   \[\text{var}(\hat{f}_h) \approx \frac{\sigma^2 V}{n h}\]

3.9 Estimación por intervalos

Los estimadores kernel pueden verse como promedios ponderados, por lo que el Teorema del Límite Central es aplicable y la distribución de \(f_\lambda\) será asintóticamente normal. Un posible enunciado para este resultado sería:

Proposición. Si los \(\epsilon_i\) son independientes e idénticamente distribuidos con \(E(\epsilon_i) = 0\) y \(\text{var}(\epsilon_i) = \sigma^2 < \infty\) y si \(n \to \infty\) y \(\lambda \to 0\) de tal modo que \(n\lambda \to \infty\), entonces para \(x \in [0, 1]\) se satisface que

\[ \frac{f_\lambda(x) - E[f_\lambda(x)]}{\sqrt{\text{var}[f_\lambda(x)]}} \xrightarrow{d} N[0, 1] \]

donde \(N[0, 1]\) representa una variable aleatoria normal estándar.

En el enunciado se acude únicamente al hecho de que \(f_\lambda\) es un promedio y el resultado no depende de suavización alguna, porque \(f_\lambda(x)-E[f_\lambda(x)]\) no depende de \(f\). El problema es más complejo cuando nos enfrentamos al estudio de la normalidad asintótica de la cantidad \(f_\lambda(x) - f(x)\), que sí requiere del uso de suavización. En este caso acudiremos al resultado:

\[ f_\lambda(x) - f(x) = \{f_\lambda(x) - E[f_\lambda(x)]\} + \{E[f_\lambda(x)] - f(x)\} \quad (3.25) \]

Si dividimos ambos lados de la expresión (3.25) entre \(\text{var}[f_\lambda(x)]\), entonces el problema de convergencia se reduce al segundo término de la derecha, que tiende a cero asintóticamente si \(n\lambda^5 \to 0\), por lo que podemos expresar el siguiente nuevo enunciado:

Proposición. Bajo las mismas condiciones del enunciado anterior, con \(f \in C^2[0, 1]\) y \(n\lambda^5 \to 0\), entonces:

\[ \frac{f_\lambda(x) - f(x)}{\sqrt{\text{var}[f_\lambda(x)]}} \xrightarrow{d} N[0, 1] \]

Llegamos entonces a un primer intervalo de confianza para \(f(x)\):

\[ f_\lambda(x) \pm Z_{\alpha/2}\hat{\sigma}\sqrt{\sum_{i=1}^n w_{in}^2(x, \lambda)} \quad (3.26) \]

Para el caso del estimador de Gasser-Müller tenemos que:

\[ w_{in}(x, \lambda) = \lambda^{-1} \int_{s_{i-1}}^{s_i} K(\lambda^{-1}(x-s))ds, \quad i = 1, \ldots, n \]

Y para el caso del estimador de Nadaraya-Watson,

\[ w_{in}(x, \lambda) = \frac{K(\lambda^{-1}(x - x_i))}{\sum_{j=1}^n K(\lambda^{-1}(x - x_j))}, \quad i = 1, \ldots, n \]

En el intervalo (3.26), \(Z_{\alpha/2}\) es el percentil \(100(1-\alpha)\) de la normal estándar y \(\hat{\sigma}\) es la raíz cuadrada de una estimación \(\hat{\sigma}^2\) de \(\sigma^2\) basada en alguno de los estimadores de la varianza en el modelo (2.1) presentados en la sección 2.2.1.

La principal limitación del intervalo (3.26) es que se construye a partir de un estimador \(f_\lambda\) de \(f\) que se sabe que no es insesgado. Vale recordar que hemos elegido \(\lambda\) de tal manera que se minimice el riesgo, que es una combinación del sesgo y de la varianza. Una opción sería intentar construir un intervalo de confianza que considere el sesgo.

Eubank (1999) propone una solución que implica estimar \(f\) con una parámetro de suavización, digamos \(\lambda_1\), estimar luego \(f''\) usando borde en los extremos. El intervalo tomaría la siguiente forma:

\[ f_{\lambda_1}(x) - \frac{\lambda^2}{2} M_2 \hat{f}_{\lambda_2}''(x) \pm Z_{\alpha/2}\hat{\sigma}\sqrt{\sum_{i=1}^n w_{in}^2(x, \lambda)} \quad (3.27) \]

En esta propuesta se debe además adecuar el cálculo de los pesos \(w_{in}(x, \lambda)\) para considerar la forma estimación de \(f\) y tiene como dificultades la elección de más de un parámetro de suavización y el uso de funciones kernel particulares en los extremos. Por otra parte, el intervalo puede verse como si se centrara en una corrección del estimador para lograr que sea insesgado.

Bowman & Azzalini (1997) prefiere utilizar “bandas de variabilidad” que se construyen usando intervalos de confianza calculados para valores de \(x\), alrededor de \(E[f_\lambda(x)]\) usando un estimador de \(\sigma^2\).

La figura 3.7 muestra los intervalos de confianza construidos usando un estimador kernel y las bandas de confianza de Bowman-Azzalini. Ambos ajustes tienen aproximadamente el mismo número de grados de libertad equivalente

# Datos: hora y CO
hora <- 0:23
CO <- c(2.24, 1.53, 1.80, 0.87, 1.21, 1.02, 1.39, 1.90, 2.49, 3.33, 2.82, 3.34,
        3.22, 2.94, 2.13, 2.88, 1.95, 1.51, 1.54, 1.60, 1.22, 1.37, 1.96, 1.66)

# Cargar la librería sm (si no está instalada, instálala con install.packages("sm"))
library(sm)
library(ggplot2)

# Establecer el marco para dos gráficos lado a lado
par(mfrow = c(1, 2))

# Estimación kernel con intervalos de confianza al 95%
#sm.regression(x = hora, y = CO, h = 3, model = "none", 
#              display = "lines", ci = TRUE,
#              lty = c(1, 2, 2), col = c("black", "red", "red"),
#              xlab = "Hora", ylab = "CO", 
#              main = "Intervalo de Confianza (95%)")

# Bandas de variabilidad al 95%
#sm.regression(x = hora, y = CO, h = 3, model = "none", 
#              display = "lines", ci = TRUE, ci.type = "variability",
#              lty = c(1, 3, 3), col = c("black", "blue", "blue"),
#              xlab = "Hora", ylab = "CO", 
#              main = "Bandas de Variabilidad (95%)")

sm.regression(x = hora, y = CO, h = 3, display = "se", col = c("black", "red"),
              xlab = "Hora", ylab = "CO")


datos <- data.frame(hora = hora, CO = CO)
ggplot(datos, aes(hora, CO)) + 
  geom_point() + 
  geom_smooth(method = "loess", se = TRUE, level = 0.95) +
  labs(title = "Regresión no paramétrica con bandas de confianza")

# Datos: hora y CO
hora <- 0:23
CO <- c(2.24, 1.53, 1.80, 0.87, 1.21, 1.02, 1.39, 1.90, 2.49, 3.33, 2.82, 3.34,
        3.22, 2.94, 2.13, 2.88, 1.95, 1.51, 1.54, 1.60, 1.22, 1.37, 1.96, 1.66)

# Cargar la librería sm (si no está instalada, instálala con install.packages("sm"))
library(sm)

# Establecer el marco para dos gráficos lado a lado
par(mfrow = c(1, 2))


sm.regression(x = hora, y = CO, h = 3, display = "se", col = c("black", "red"),
              xlab = "Hora", ylab = "CO")

4 Estimación spline

El término spline proviene del nombre de una herramienta utilizada por ingenieros y arquitectos navales para trazar curvas suaves sobre planos de navios, bajo este concepto de crear curvas suaves que pasan a través de un conjunto de puntos algunos matemáticos comenzaron a trabajar en ello, para poder representar curvas de manera más flexible y precisa que con las técnicas anteriores, como las interpolaciones lineales..

4.0.1 Interpolación y suavización spline

La interpolación es una técnica matemática que nos permite estimar valores desconocidos en función de datos conocidos. Básicamente, consiste en construir una función o curva que pase a través de puntos específicos para luego calcular valores intermedios con base en esa función.

Existe un diferente metodos de interpolación, pero para nuestro estudio hablaremos de la interpolación polinomica, en donde, si se dispone de un conjunto \(A\) de \(k\) puntos en \(\mathbb{R}^2\) entonces se podrá ajustar a los puntos exactamente un polinomio de grado \(k-1\). Es decir, con dos puntos se puede ajustar una recta y con cuatro puntos una cúbica. En esencia si se cuenta con \(K\) puntos entonces el polonomio que se quiere interpolar tendra la forma

\[\pi_{k-1}=a_0+a_1x+a_3x^2+\dots+a_kx^{k-1}.\] Para determinar completamente el polinomio \(\pi_{k-1}\) basta con encontrar los coeficientes \(a_0,a_1,\dots,a_k\) para lo cual se utilizan las coordenadas \((x_i,y_i),i=1,2,\dots,k\) de los k puntos de A. Para esto se requieren \(k\) ecuaciónes, una para cada coeficiente, las que se construyen usando los \(k\) puntos de \(A\).

Por ejemplo, supongamos que tenemos

\[A=\{(1,2);(2,5);(3,4);(4,7)\}\] Y sea \(\pi_3=f(x)=a_0+a_1x+a_2x^2\) el polinomio cuya gráfica pasa por estos cuatro puntos. Para identificar el polinomio \(\pi_3\) usando interpolación, basta con encontrar los valores de las constantes \(a_0,a_1,a_2\).

Al resolver obtenemos que los valores son: \((-13,23.66,-10)\) y el gráfico del polinomio aparece a continuación.

# Datos
x <- c(1, 2, 3, 4)
y <- c(2, 5, 4, 7)

# Ajuste del polinomio de grado 3
modelo <- lm(y ~ poly(x, degree = 3, raw = TRUE))

# Ver coeficientes del polinomio
coeficientes <- c(modelo$coefficients[1],modelo$coefficients[2],modelo$coefficients[3])
coeficientes

##                      (Intercept) poly(x, degree = 3, raw = TRUE)1 
##                        -13.00000                         23.66667 
## poly(x, degree = 3, raw = TRUE)2 
##                        -10.00000

# Secuencia de puntos para graficar el polinomio
x_vals <- seq(1, 4, by = 0.01)
y_vals <- predict(modelo, newdata = data.frame(x = x_vals))

# Graficar puntos originales y curva interpolada
plot(x, y, col = "red", pch = 19, main = "Interpolación Polinómica de Grado 3", xlab = "x", ylab = "y")
lines(x_vals, y_vals, col = "blue", lwd = 2)

Una propiedad llamativa de un polinomio cúbico es que tiene un único punto de inflexión. Si se desea construir una curva que tiene más puntos de inflexión, una solución posible es resolver un sistema más grande de ecuaciones lineales. Otra solución posible en este caso sería construir la interpolación uniendo varios polinomios cúbicos.

Para lograrlo se escogen varios puntos llamados comúnmente . Sean \((p_1, p_2, \ldots, p_n)\) los puntos de control y \((f_1, f_2, \ldots, f_{n-1})\) las cúbicas que se usarán para interpolar \(f\).

El propósito es ajustar una cúbica entre cada par de puntos de control \((p_i, p_{i+1})\). Si se denota \(f_i(x) = a_{i0} + a_{i1}x + a_{i2}x^2 + a_{i3}x^3\) a la cúbica que une estos dos puntos de control.

Para construir estas ecuaciones es necesario imponer algunas condiciones mínimas a las cúbicas \(f_i\) que se unirán entre sí para construir el polinomio \(f\).

4.0.1.1 Estimación

Un problema igualmente interesante desde el punto de vista numérico pero más interesante desde el punto de vista estadístico es aquel en el cual se sabe que los datos disponibles son mediciones sujetas a error y que aunque se cree que la relación entre \(X\) y \(Y\) es funcional, no necesariamente los puntos observados están sobre la curva que relaciona las dos variables, como se ve en la siguiente grafica

Entonces, el propósito de la estimación es proponer una función \(\hat{f}\) que permita aproximar \(f\).

En realidad en ambos casos, en la interpolación y en la estimación, se obtienen curvas aproximadas. La diferencia está en que, por una parte, no es posible utilizar interpolación con los puntos y que a través de los métodos de estimación es posible asociar a las curvas construidas una medida probabilística de qué tan buena será la aproximación alcanzada.

Una posible solución estadística es la regresión polinómica de la sección 2.3.1, en la que se asume que la función subyacente \(f\) es una función continua que se aproxima mediante un polinomio de Taylor y que las \(n\) observaciones disponibles provienen de tal función \(f\), de tal manera que los valores \(y_i\) se separan de los de \(f(x_i)\) por una cantidad aleatoria \(\varepsilon_i\).

4.0.1.2 Estimación splaine

Una función spline es una función suave a tramos que se utiliza para aproximar a a un conjunto de datos atraves de splines, estos unen varios polinomios de bajo grado en distintos intervalos, asegurando que la transición entre ellos sea suave.

Supongamos que tenemos un conjunto de números reales \(x_1, \ldots, x_n\) en un intervalo \([a, b]\), tales que \[ a < x_1 < x_2 < \ldots < x_n < b. \]

Una función \(s\) definida en \([a, b]\) es un spline cúbico si cumple las siguientes dos condiciones:

1. \(s\) es una función cúbica en cada uno de los intervalos \((a, x_1), (x_1, x_2), \ldots, (x_n, b)\).
2. Las cúbicas se unen en los puntos \(x_i\) de tal manera que \(s\) y sus dos primeras derivadas son continuas en cada \(x_i\), y por lo tanto en todo el intervalo \([a, b]\).

A los puntos \(x_i\) los llamaremos nodos.

Definición Un spline cúbico en \([a, b]\) se llama un spline cúbico natural si se cumple que las dos primeras derivadas de \(s\) son iguales a cero en los puntos extremos \(a\) y \(b\). Estas condiciones se conocen como condiciones de acotamiento natural, e implican que \(s\) es lineal en los dos intervalos extremos \((a, x_1)\) y \((x_n, b)\).

En la estimación usando series de Fourier acordamos que la función de regresión está en \(L_2[0, 1]\). Añadiremos ahora la condición de que las dos primeras derivadas de \(f\), que denotaremos \(f'\) y \(f''\), también pertenezcan a \(L_2[0, 1]\). A este espacio de funciones lo llamaremos el espacio de Sobolev de orden 2 y lo denotaremos por \[ W_2^2[0, 1]. \]

Una medida natural de suavidad asociada con una función \(f \in W_2^2[0, 1]\) es \[ \int_0^1 \left( f''(x) \right)^2 \, dx, \] mientras que una medida de bondad de ajuste de los datos al modelo es la suma de cuadrados del error: \[ \frac{1}{n} \sum_{i=1}^{n} \left( y_i - f(x_i) \right)^2. \] Esto implica que una medida de la calidad de un estimador de \(f\) podría basarse en la siguiente suma convexa: \[ (1 - q)\frac{1}{n} \sum_{i=1}^{n} \left( y_i - f(x_i) \right)^2 + q \int_0^1 \left( f''(x) \right)^2 dx, \] con \(0 < q < 1\).

Si hacemos \(\lambda = \dfrac{q}{1 - q}\), la elección del estimador de \(f\) es equivalente a elegir \(f_\lambda\) que minimice la siguiente expresión:

\[\begin{equation} \frac{1}{n} \sum_{i=1}^{n} \left( y_i - f(x_i) \right)^2 + \lambda \int_0^1 \left( f''(x) \right)^2 dx, \tag{4.1} \end{equation}\] sobre todas las funciones \(f \in W_2^2[0, 1]\). A este estimador \(f_\lambda\) lo llamaremos un estimador spline de \(f\).

De la expresión (4.1) se sigue que:

• Si \(\lambda\) es muy grande, la estimación de la función de regresión será super-suavizada.
• Si \(\lambda\) es muy pequeño, se obtiene un estimador que interpola exactamente los datos.

Eubank en su libro “Nonparametric Regression and Spline Smoothing” demostro que la solución a este problema de optimización es única, y corresponde al estimador: \[ f_\lambda = \sum_{i=1}^{n} \beta_{\lambda i} f_i, \] donde los \(\beta_{\lambda i}\) son coeficientes reales, y \(\{f_i\}\) es un conjunto de funciones base..

Donde \(\boldsymbol{\beta}_\lambda = (\beta_{\lambda 1}, \beta_{\lambda 2}, \ldots, \beta_{\lambda n})^\top\) es la única solución con respecto a \[ \mathbf{c} = (c_1, c_2, \ldots, c_n)^\top \] del sistema de ecuaciones: \[\begin{equation} (X^\top X + n\lambda \, \Omega) \, \mathbf{c} = X^\top \mathbf{y} \end{equation}\] donde: \[ X = \{f_j(x_i)\}_{i,j=1,2,\ldots,n}, \quad \mathbf{y} = (y_1, y_2, \ldots, y_n)^\top, \] y \[ \Omega = \left[ \int_0^1 f_i''(x) f_j''(x) \, dx \right]_{i,j=1,2,\ldots,n}. \]

Las funciones \(\{f_j\}_{j=1}^{n}\) forman una base del conjunto de splines cúbicos naturales.

Si usamos la base de splines cúbicos naturales, entonces se tendría que el vector de valores estimados es \[ \mathbf{f}_\lambda = \begin{pmatrix} f_\lambda(x_1), f_\lambda(x_2), \ldots, f_\lambda(x_n) \end{pmatrix}^T = S_\lambda \mathbf{y}, \] donde \[\begin{equation} S_\lambda = X (X^T X + n\lambda \Omega)^{-1} X^T. \tag{4.2} \end{equation}\]

La elección del parámetro de suavización \(\lambda\) se hace usualmente con el estimador de validación cruzada generalizada GCV, usando la matriz \(S_\lambda\).

En R existe la función smooth.spline del paquete stats que selecciona \(\lambda\) y realiza la estimación.

#datos
hora_reescalada <- (0:23)/23
co <- c(2.24, 1.53, 1.80, 0.87, 1.21, 1.02, 1.39, 1.90, 2.49, 3.33, 
        2.82, 3.34, 3.22, 2.94, 2.13, 2.88, 1.95, 1.51, 1.54, 1.60, 
        1.22, 1.37, 1.96, 1.66)

# Ajuste sin spar
spline <- smooth.spline(hora_reescalada, co)
print(spline)

## Call:
## smooth.spline(x = hora_reescalada, y = co)
## 
## Smoothing Parameter  spar= 0.5058871  lambda= 0.0002605305 (14 iterations)
## Equivalent Degrees of Freedom (Df): 7.306994
## Penalized Criterion (RSS): 1.577585
## GCV: 0.1358735

# Gráfico
plot(hora_reescalada, co, col="blue", pch=19, 
     main="Ajuste con smooth.spline")
lines(predict(spline, seq(0,1,0.01)), col="red", lwd=2)

Aunque tambien la función te permite modificar el parámetro spar en la función smooth.spline() de R controla el grado de suavizado del ajuste del spline cúbico suavizado, este valor esta entre 0 y 1, donde 1 es la maxima suavización.

#Datos
hora_reescalada <- (0:23)/23
co <- c(2.24, 1.53, 1.80, 0.87, 1.21, 1.02, 1.39, 1.90, 2.49, 3.33, 
        2.82, 3.34, 3.22, 2.94, 2.13, 2.88, 1.95, 1.51, 1.54, 1.60, 
        1.22, 1.37, 1.96, 1.66)

# Ajuste sin spar
spline_1 <- smooth.spline(hora_reescalada, co,spar = 0.95)
spline_2 <- smooth.spline(hora_reescalada, co,spar = 0.4)

# Gráfico
plot(hora_reescalada, co, col="blue", pch=19)
lines(predict(spline_1, seq(0,1,0.01)), col="red", lwd=2)
lines(predict(spline_2, seq(0,1,0.01)), col="purple", lwd=2)
legend("topright", legend = c("spar = 0.95", "spar = 0.4"),
       col = c("red", "purple"), lwd = 1, bty = "n", cex = 1)

4.0.2 Estimación spline por minimos cuadrados

Las propiedades asintóticas de los estimadores spline que hemos discutido hasta aquí se derivan del hecho de que estos estimadores pueden representarse como estimadores kernel, de tal manera que sus propiedades asintóticas son similares.

Existe, sin embargo, otro posible estimador spline que no puede representarse como un estimador kernel. Esta nueva forma de estimación se apoya en el modelo de regresión polinómica discutido. En ese caso, representamos el modelo (2.1) de la siguiente manera:

\[\begin{equation} f(x_i) = \sum_{j=0}^{\lambda} \theta_j f_j(x_i) + r_\lambda(x_i), \quad i = 1, \ldots, n \tag{4.3} \end{equation}\]

donde \(\lambda\) cumple el papel de parámetro de suavización, y consiste en el número de términos que se conservan en la sumatoria.

Si modificamos la notación y llamamos \(m\) a lo que en la expresión (4.3) hemos denotado como \(\lambda\), una de las representaciones posibles para \(r_m(x)\) sería:

\[\begin{equation} r_m(x) = \frac{1}{m - 1} \int_0^1 f^{(m)}(\xi) (x - \xi)^{m-1}_+ \, d\xi \end{equation}\]

Si \(r(x_1), \ldots, r(x_n)\) son muy pequeños, entonces el modelo de regresión polinómica sería adecuado para un conjunto particular de datos que tenga estas características. Pero si este no es el caso, una posible generalización es la , que se apoya en la siguiente aproximación a la integral:

\[\begin{equation} r_m(x) \approx \sum_{j=1}^{k} \delta_j (x - \xi_j)^{m-1}_+ \end{equation}\]

para un conjunto de constantes \(\delta_1, \ldots, \delta_k\) y un conjunto de puntos \(0 < \xi_1 < \ldots < \xi_k < 1\), de tal manera que una nueva aproximación a la función de regresión lucirá ahora así:

\[\begin{equation} f(x) = \sum_{j=1}^{\lambda} \theta_j x^{j-1} + \sum_{j=1}^{k} \delta_j (x - \xi_j)^{m-1}_+ \end{equation}\]

Sea \(\{\xi_1, \ldots, \xi_k\}\) un conjunto dado de puntos (nodos). Una posible solución para estimar \(f\) sería estimar por mínimos cuadrados los \(m\) coeficientes \(\{\theta_j\}_{j=1}^{m}\) y los \(k\) coeficientes \(\{\delta_j\}_{j=1}^{k}\).

Sea entonces el vector de parámetros:

\[ \boldsymbol{\beta} = (\theta_1, \ldots, \theta_m, \delta_1, \ldots, \delta_k)^\top \]

El estimador spline por mínimos cuadrados de \(f\) será:

\[\begin{equation} f_\lambda(x) = \sum_{j=1}^{m+k} \beta_{\lambda j} f_j(x) \end{equation}\]

Un problema aún en curso es la elección de \(m\) y de \(k\). A los \(k\) valores de \(\xi\) se les llama nodos, y la solución más sencilla es hacer uso de inspección visual de los datos, cuando sea posible.

Si \(f\) cambia muy rápidamente, entonces deberemos usar más nodos. La elección de \(m\) depende de si se usará un spline lineal, cuadrático, cúbico, etc., y los \(k\) nodos se ubicarán consecuentemente.

En R existe la funcíon bs para ajustar splines con mínimos cuadrados ordinarios, la función genera una base de B-splines (Basis Splines) para ajustar modelos de splines.

library(splines)
# Datos
x <- (0:23)/23  
y <- c(2.24, 1.53, 1.80, 0.87, 1.21, 1.02, 1.39, 1.90, 
       2.49, 3.33, 2.82, 3.34, 3.22, 2.94, 2.13, 2.88, 
       1.95, 1.51, 1.54, 1.60, 1.22, 1.37, 1.96, 1.66)

# Definir los nodos 
knots <- quantile(x, probs = seq(0.1, 0.9, length.out = 5))

# Crear matriz  B-splines cúbicos
X <- bs(x, knots = knots, degree = 3, intercept = TRUE)

# Ajustar modelo por mínimos cuadrados
fit <- lm(y ~ X - 1)  # EL  intercepto ya esta  incluido en bs()

# Crear puntos para la curva ajustada
x_new <- seq(0, 1, length.out = 100)
X_new <- predict(bs(x, knots = knots, degree = 3, intercept = TRUE), x_new)
y_pred <- X_new %*% coef(fit)

# Graficar resultados
plot(x, y, pch = 19, col = "blue", main = "B-splines")
lines(x_new, y_pred, col = "red", lwd = 2)
points(knots, predict(fit, newdata = list(X = predict(bs(x, knots = knots, degree = 3, intercept = TRUE), knots))), 
       pch = "x", col = "green", cex = 1.5)
legend("topright", legend = c("Ajuste", "Nodos"), 
       col = c("red", "green"), pch = c(NA, "x"), lwd = c(2, NA))

5 Modelos aditivos Generealizados

Los Modelos Aditivos Generalizados (GAMs) extienden los modelos lineales permitiendo relaciones no lineales entre predictores y la variable respuesta, manteniendo una estructura aditiva.

Para una variable respuesta \(Y\) con distribución en la familia exponencial y predictores \(X_1,\ldots,X_p\), el GAM se define como:

\[ g\left(\mathbb{E}[Y|X]\right) = \beta_0 + \sum_{j=1}^p f_j(X_j) \]

donde:

• \(g:\mathbb{R}\to\mathbb{R}\) es la función de enlace
• \(f_j\) son funciones suaves (típicamente splines cúbicos)
• \(Y \sim F(\mu,\phi)\) con \(\mu = \mathbb{E}[Y|X]\)

\[ \begin{aligned} \text{Función de enlace:} \quad & \eta = g(\mu) \\ \text{Función de varianza:} \quad & \text{Var}(Y) = V(\mu)\phi \\ \text{Penalización:} \quad & \int \left[f_j''(x)\right]^2 dx \end{aligned} \]

Se maximiza la log-verosimilitud penalizada, para mas detalles José (2022)

\[ \ell_p(\boldsymbol{\theta}) = \ell(\boldsymbol{\theta}) - \frac{1}{2}\sum_{j=1}^p \lambda_j \int \left[f_j''(x)\right]^2 dx \]

donde \(\boldsymbol{\theta}\) incluye \(\beta_0\) y los parámetros de las funciones suavizadas \(f_j\).

Familias Comunes

Familia	Función de enlace	Aplicación
Gaussiana	Identidad	Datos continuos
Gamma	Log	Datos positivos
Poisson	Log	Conteos
Binomial	Logit	Proporciones

A modo de compración:

Característica	GLM	GAM
Forma funcional	Lineal (\(\beta_j X_j\))	No lineal (\(f_j(X_j)\))
Término error	\(\varepsilon\) explícito	Implícito en distribución

Para los datos en la sección @(Modelos paramétricos)

library(ggplot2)
library(mgcv)

set.seed(123)
x <- seq(1, 10, length.out = 100)
y <- exp(0.3 * x + rnorm(100, sd = 0.2))
datos <- data.frame(x = x, y = y)

modelo1 <- gam(y ~ s(x), data = datos)
modelo2 <- gam(y ~ s(x, k = 15), family = Gamma(link = "log"), method = "REML", data = datos)


AIC(modelo1, modelo2)

##               df      AIC
## modelo1 6.022066 404.9719
## modelo2 3.003534 282.5270

summary(modelo1)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## y ~ s(x)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   7.3333     0.1771   41.41   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##        edf Ref.df     F p-value    
## s(x) 4.022  4.965 208.1  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.913   Deviance explained = 91.6%
## GCV = 3.3017  Scale est. = 3.1359    n = 100

summary(modelo2)

## 
## Family: Gamma 
## Link function: log 
## 
## Formula:
## y ~ s(x, k = 15)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.68453    0.01851   91.01   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##        edf Ref.df    F p-value    
## s(x) 1.002  1.004 1868  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.912   Deviance explained = 94.9%
## -REML = 144.44  Scale est. = 0.034258  n = 100

par(mfrow = c(3, 2))
gam.check(modelo1, old.style = TRUE)

## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 6 iterations.
## The RMS GCV score gradient at convergence was 3.073791e-07 .
## The Hessian was positive definite.
## Model rank =  10 / 10 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##        k'  edf k-index p-value
## s(x) 9.00 4.02    1.03    0.68

par(mfrow = c(2, 2))
gam.check(modelo2, old.style = TRUE)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 9 iterations.
## Gradient range [-7.231766e-05,0.0007281679]
## (score 144.4377 & scale 0.03425775).
## Hessian positive definite, eigenvalue range [7.344579e-05,49.55541].
## Model rank =  15 / 15 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##      k' edf k-index p-value
## s(x) 14   1    1.03    0.59

par(mfrow = c(1, 2)) 
plot(modelo1, shade = TRUE, main = "Modelo 1: Básico")
plot(modelo2, shade = TRUE, main = "Modelo 2: Gamma-log")

6 Apendice

6.0.1 Espacio \(L_2[0,1]\)

Definición: EL espacio \(L_2[0,1]\) se define como el espacio funcional de todas las funciones cuadrado integrables en [0, 1].

Sean \(f_1\) y \(f_2\) dos funciones en \(L_2[0, 1]\). Entonces, la norma de \(f_i\) se define como:

\[ \|f_i\| = \left( \int_0^1 f_i^2(x) \, dx \right)^{1/2}, \quad i = 1, 2 \] mientras que el producto interno de \(f_1\) y \(f_2\) es:

\[ <f_1, f_2 > = \int_0^1 f_1(x) f_2(x) \, dx \]

6.0.1.1 Sucesiones ortonormales completas (CONS)

En un espacio de dimensión finita, por ejemplo en \(\mathbb{R}^p\), siempre es posible representar cualquier elemento del espacio usando una combinación lineal de elementos en una base del espacio.

Pero el espacio \(L_2[0, 1]\) es un espacio de dimensión infinita, por lo que las siguientes definiciones se requieren para poder concluir que, en efecto, es posible representar cualquier elemento de \(L_2[0, 1]\) a partir de un conjunto de elementos que formen una base de \(L_2[0, 1]\).

1. Dos funciones \(f_1, f_2 \in L_2[0, 1]\) se dicen ortogonales si \(\langle f_1, f_2 \rangle = 0\).
   La ortogonalidad de \(f_1\) y \(f_2\) se denota \(f_1 \perp f_2\).

2. Una sucesión de funciones \(\{f_j\}_{j=1}^{\infty} \subset L_2[0, 1]\) se dice ortonormal
   si las \(f_j\) son ortogonales por parejas y \(\|f_j\| = 1\) para todo \(j\).

3. Una sucesión de funciones \(\{f_j\}_{j=1}^{\infty} \subset L_2[0, 1]\) se dice que es una
   sucesión ortonormal completa (CONS, por sus siglas en inglés) si
   \(f \perp f_j\) para todo \(j\) implica que \(f = 0\).

Las definiciones de ortogonalidad y ortonormalidad en \(L_2[0, 1]\) son similares a las mismas definiciones en \(\mathbb{R}^p\), mientras que la definición de una CONS implica que la única función ortogonal a todas las funciones \(f_j\) es la función cero.

6.0.1.2 Coeficientes de Fourier

La siguinete proposición nos permite concluir que, cualquier CONS provee una base de \(L_2[0, 1]\). Por lo que una función \(f\) de \(L_2[0, 1]\) puede representarse como una combinación lineal de una colección de funciones \(\{f_j\}_{j=1}^\infty\).

Proposicion Sea \(\{f_j\}_{j=1}^{\infty}\) una CONS de \(L_2[0,1]\) y sea \(f\) una función de \(L_2[0,1]\). Si definimos \[\begin{equation} f_\lambda(x)=\sum_{i=1}^nK(x, x_i; \lambda)y_i \tag{6.1} \end{equation}\] entonces \(\sum_{j=1}^{\lambda} \beta_j f_j\) es la mejor aproximación a \(f\) en el sentido de que \[ \left\| f - \sum_{j=1}^{\lambda} \beta_j f_j \right\| \leq \left\| f - \sum_{j=1}^{\lambda} b_j f_j \right\| \] para todo \(b = (b_1, b_2, \dots, b_\lambda)^T \in \mathbb{R}^\lambda\). Más aún, cuando \(\lambda \to \infty\), \[ \left\| f - \sum_{j=1}^{\lambda} \beta_j f_j \right\| \to 0. \]

Llamaremos coeficientes generalizados de Fourier a los coeficientes \(\beta_j\) de la expresión (6.1).

Finalmente, se sabe que los coeficientes generalizados de Fourier satisfacen la relación de Parseval. \[\sum_{j=1}^\infty \beta_j^2=\|f\|.\]

6.0.1.3 Diferencias de orden

Definicón Dada una secuencia de datos \(\{y_i\}_{i=1}^n\) una diferencia de orden \(r\) con pesos \(\{d_k\}_{k=-m_1}^{m_2}\), se calcula como: \[ \Delta_l=\sum_{k = -m_1}^{m_2} d_k y_{k + l}. \] donde los \(d_k\) satisfacen: \(\sum_{k=-m_1}^{m_2} d_k = 0\) (elimina tenedencias de bajo grado) y \(\sum_{k=-m_1}^{m_2} d_k^2 = 1\) (mantiene la escala de la varianza).

6.0.2 Espline Cubico

Dado un conjunto de nodos: \[ a = x_0 < x_1 < x_2 < \cdots < x_n = b, \] un spline cúbico \(S(x)\) está definido como: \[ S(x) = \begin{cases} P_1(x), & x \in [x_0, x_1] \\ P_2(x), & x \in [x_1, x_2] \\ \vdots \\ P_n(x), & x \in [x_{n-1}, x_n] \end{cases} \] donde cada \(P_i(x)\) es un polinomio cúbico de la forma: \[ P_i(x) = a_i + b_i(x - x_{i-1}) + c_i(x - x_{i-1})^2 + d_i(x - x_{i-1})^3 \]

Referencias

Benedetti, Jacqueline K. 1975. “KERNEL ESTIMATION OF REGRESSION FUNCTIONS.”

Berrendero, José R. 2023. “Estimadores Núcleo En Regresión No Paramétrica.” 2023. https://verso.mat.uam.es/~joser.berrendero/libro-est/estimadores_nucleo.html.

Cleveland, William S. 1979. “Robust Locally Weighted Regression and Smoothing Scatterplots.” Journal of the American Statistical Association 74 (368): 829–36.

Eubank, Randall L. 1999. Nonparametric Regression and Spline Smoothing. CRC press.

Green, Peter J, and Bernard W Silverman. 1993. Nonparametric Regression and Generalized Linear Models: A Roughness Penalty Approach. Crc Press.

José, Clemente García Francisco. 2022. “MODELOS LINEALES GENERALIZADOS y ADITIVOS GENERALIZADOS.” Master’s thesis, Universidad de Extremadura. https://dehesa.unex.es/bitstream/10662/16739/1/TFGUEX_2022_Clemente_Garcia.pdf.

Nadaraya, Elizbar A. 1964. “On Estimating Regression.” Theory of Probability & Its Applications 9 (1): 141–42.

Ochoa, Javier Olaya. 2020b. “MÉTODOS DE REGRESIÓN NO PARAMÉTRICA.” Prospectiva. Escuela de Estadística. https://doi.org/10.25100/peu.505.

———. 2020a. “MÉTODOS DE REGRESIÓN NO PARAMÉTRICA.” Prospectiva. Escuela de Estadística. https://doi.org/10.25100/peu.505.

Priestley, Maurice B, and MT Chao. 1972. “Non-Parametric Function Fitting.” Journal of the Royal Statistical Society Series B: Statistical Methodology 34 (3): 385–92.

Wahba, Grace. 1990. Spline Models for Observational Data. SIAM.

Watson, Geoffrey S. 1964. “Smooth Regression Analysis.” Sankhyā: The Indian Journal of Statistics, Series A, 359–72.