Elección de lambda, regresión no paramétrica
Julián David Téllez
9/4/2020
Estimación de \(\sigma^2\) y elección del \(\lambda\) óptimo
Seleccionar uno de los cinco días del ejercicio anterior.
Estimar la \(\sigma^2\) del modelo usando el estimador de \(Rice\).
Usar el estimador de \(Rice\) en el estimador \(UBRE\) para elegir el mejor \(\lambda\).
Calcular los estimadores de validación cruzada (\(CV\)) y de validación cruzada generalizada (\(GCV\)) para elegir el mejor \(\lambda\).
Librerias necesarias
library(readxl)
library(dplyr)
library(lubridate)
library(ggplot2)
library(gridExtra)
library(purrr)Carga de Datos
Para efectos de este ejercicio Se hace uso del \(dataset\) \(DatosLaFlora\) el cual contiene información sobre mediciones horarias de \(NO_2\) y \(O_3\) en (\(\mu g/m3\)) y sobre algunas variables de clima en el norte de Cali (Colombia). Se dispone además de un campo que informa sobre el día y la hora a la que se tomó la medición.
Datos <- read_excel(ruta)Ajuste de los datos
Lo que se busca en esta sección es conformar un \(dataframe\) con las variables (\(NO_2\), \(dia\) y \(hora\_normada\)) del dataset \(DatosLaFlora\), las variables \(dia\) y \(hora\_normada\) se obtienen a partir de la variable \(Fecha \& Hora\); \(hora\_normada\) se usa para restringirnos al intervalo \([0,1]\), este espacio de las funciones cuadrado integrables se denota en general como \(L_2[0,1]\) el cual satisface las propiedades generales de un espacio vectorial.
La transformación se hace como sigue:
para cada hora \(i \in [1,24]\) se aplica la siguiente ecuación para normalizar:
\[x_i=\frac{(2_i-1)}{2n}, \hspace{1 cm} \forall i=1,2,...,24\]
x <- Datos %>%
dplyr::select(`Fecha & Hora`, NO2) %>%
mutate(dia = day(`Fecha & Hora`)) %>%
mutate(fecha = as.Date(Datos$`Fecha & Hora`)) %>%
filter(dia %in% 10) %>%
mutate(hora = row_number()) %>%
mutate(hora_normada = (2*hora - 1)/(2*24)) %>%
dplyr::select(fecha, hora_normada, NO2)Obtenemos entonces el siguiente \(tibble\):
## # A tibble: 24 x 3
## fecha hora_normada NO2
## <date> <dbl> <dbl>
## 1 2011-08-10 0.0208 63.8
## 2 2011-08-10 0.0625 46.1
## 3 2011-08-10 0.104 47.7
## 4 2011-08-10 0.146 39.9
## 5 2011-08-10 0.188 35.2
## 6 2011-08-10 0.229 31.2
## 7 2011-08-10 0.271 26.8
## 8 2011-08-10 0.312 22.9
## 9 2011-08-10 0.354 37.0
## 10 2011-08-10 0.396 35.3
## # ... with 14 more rowsExploración de los datos
Generamos un gráfico de puntos con las mediciones de \(NO_2\) vs \(hora\_normada\) para observar el comportamiento de los datos a lo largo del día \(10\).
ggplot()+
geom_point(data = x, aes(x = hora_normada, y = NO2)) +
labs(y = expression(NO[2]))
Ajuste del modelo
Hacemos uso del modelo del ejercicio anterior donde usamos el estimador de series de Fourier y una base \((CONS)\) de cosenos:
\[\mu_{\lambda}(x)=\sum_{j=1}^\lambda \beta_{\lambda j}\sqrt{2}cos[(j-1)\pi x]\] ¿Cuál valor de \(\lambda\) seria entonces “la mejor” elección?
Estimacion de \(\hat{\sigma}^2\)
Para la estimación de \(\hat{\sigma}^2\) usaremos un estimador de \(\sigma^2\) que no depende de la elección de \(\lambda\), la idea original fue propuesta por John Rice en 1984, y lo llamaremos \(estimador \hspace{0.1 cm} Rice\) el cual se define como sigue:
\[\sigma_{R}^2=\frac{1}{2(n-1)}\sum_{i=2}^n(y_i-y_{i-1})^2\] \(\sigma_{R}^2\) depende en esencia de las contribuciones de los errores aleatorios que contienen toda la información acerca de \(\sigma^2\).
n = nrow(x)
y = pull(x, NO2)
sigma.rice <- 1/(2*(n-1))*sum((y - lag(y, k=1))^2, na.rm = T)Tenemos entonces que \(\sigma_{R}^2=\)
## [1] 64.11709Selección del \(\lambda\)
En esta sección lo que buscamos es encontrar el \(\lambda\) que minimice el riesgo \(R(\lambda)\). Utilizaremos el estimador \(UBRE\) para identificar el valor de \(\lambda\) más adecuado para la estimación de \(\mu\) mediante series de cosenos.
El estimador \(UBRE\) (estimador insesgado del riesgo) con un estimador de la varianza que no dependa de la elección de \(\lambda\) seria:
\[\hat{R}(\lambda)=\frac{1}{n}RSS(\lambda)+\frac{2}{n}\hat{\sigma}^2tr[\textbf{S}_{\lambda}]-\hat{\sigma}^2\]
Donde \(\lambda \in \{1,2,\dots,24\}\), es el número de funciones \(f_i\) que conservaremos en la estimación de \(\mu\).
Calculamos entonces para cada valor posible de \(\lambda\) el estimador \(UBRE\) del riesgo \(\hat{R}(\lambda)\). Este estimador tiene como principal limitación su dependencia del conocimiento de la varianza \(\sigma^2\), por esto usaremos métodos de selección de \(\lambda\) que no dependan de \(\sigma^2\) como son los propuestos por Green y Silverman (2000), los cuales se definen acontinuacion.
Validación cruzada \(CV_{GS}\)
\[CV_{GS}(\lambda)=\frac{1}{n}\sum_{i=1}^n \left(\frac{y_i-\mu_{\lambda i}}{1-\textbf{S}_{\lambda ii}}\right)^2\] donde \(\textbf{S}_{\lambda ii}\) es el elemento \(i\) de la matriz \(\textbf{S}_{\lambda}\).
Validación cruzada generalizada \(GCV_{GS}\)
\[GCV_{GS}(\lambda)=n^{-1} \frac{\sum_{i=1}^n (y_i - \mu_{\lambda i})^2}{(1-n^{-1}tr[\textbf{S}_{\lambda}])^2}\]
Deseamos entonces construir un \(dataframe\) tomando como variable respuesta \(NO_2\) y como variables predictoras \(f_1, f_2, \dots, f_j\) donde \(j=\lambda\) y \(f\) es la base de cosenos \((CONS)\) que elegimos previamente.
Una vez construido el \(data.frame\) procedemos con la estimación por mínimos cuadrados de los coeficientes de Fourier haciendo uso de la función lm(). Esta estimación debe realizarse para cada valor de \(\lambda\).
Con el resultado anterior procedemos con el cálculo de \(\textbf{S}_{\lambda}, tr[\textbf{S}_{\lambda}], UBRE, CV\) y \(GCV\).
base_cons <- function(x,j){
sqrt(2)*cos((j-1)*pi*x)
}
#-------------------------------
lambda.select <- function(x, lambda, salida=1){
df <- dplyr::select(x, NO2)
hora_normada <- x$hora_normada; i <- list(); y <- list()
for (i in 2:lambda) {
y[[i]] <- base_cons(hora_normada, i)
}
y <- data.frame(matrix(unlist(y), ncol = lambda-1))
df <- bind_cols(df, y)
#-------------------------------
f.i <- df %>%
dplyr::select(contains("x")) %>%
colnames()
lambda <- lambda %>%
tibble(lambda = . )
f.i_sum <- paste(f.i, collapse = "+")
mdl_formula <- as.formula(paste("NO2", f.i_sum, sep = "~"))
mdl <- lm(mdl_formula, data = df)
#-------------------------------
fitted <- predict(mdl) %>%
tibble(fitted = . )
#-------------------------------
S = lm.influence(mdl)$hat
tr = sum(S)
#-------------------------------
UBRE <- (1/n) * sum(resid(mdl)^2) + (2/n) * sigma.rice*tr - sigma.rice %>%
tibble(ubre = .)
CV <- (1/n) * sum(((residuals(mdl) / (1-S))^2)) %>%
tibble(cv = .)
GCV <- (1/n) * ( sum(resid(mdl)^2) / (1 - (1/n) * tr)^2 ) %>%
tibble(GCV = .)
#-------------------------------
R <- bind_cols(UBRE, CV, GCV, lambda)
#-------------------------------
if(salida == 1) {
return(R)
} else {
return(fitted)
}
}
#-------------------------------
all.R <- function(x, lambda){
R <- list()
for(i in 2:lambda){
R[[i]] <- lambda.select(x,i)
}
R <- as.data.frame(t(matrix(unlist(R), ncol = lambda-1)))
names(R) <- c("UBRE","CV","GCV","LAMBDA")
return(R)
}
#-------------------------------Obtenemos finalmente un \(dataframe\) con el valor del \(\hat{R}(\lambda)\) para cada \(\lambda\) y para cada método \((UBRE, CV, GCV)\)
lambda <- 23
all.R <- all.R(x, lambda)## UBRE CV GCV LAMBDA
## 1 160.807607 260.68298 254.96154 2
## 2 139.733480 255.01389 245.31761 3
## 3 130.974007 261.69432 250.15498 4
## 4 9.184820 76.50287 74.33185 5
## 5 14.248850 83.36310 82.32426 6
## 6 19.590330 96.08913 92.29111 7
## 7 11.992016 84.01541 75.06985 8
## 8 6.169403 60.45775 56.82861 9
## 9 9.979869 60.62782 60.73289 10
## 10 15.147132 70.76860 69.83650 11
## 11 19.158601 89.08745 76.63440 12
## 12 20.730004 87.76998 73.24679 13
## 13 24.640457 86.95268 80.37662 14
## 14 29.612089 114.39383 96.58891 15
## 15 34.948796 166.30381 122.18789 16
## 16 40.124170 254.35059 157.62082 17
## 17 38.068071 89.31907 96.15241 18
## 18 42.943619 113.39058 127.68727 19
## 19 47.835601 296.26618 183.27145 20
## 20 52.509505 980.57680 282.98796 21
## 21 54.877894 1271.29394 208.36572 22
## 22 59.241387 3230.57342 269.21445 23Visualización
Finalmente con el nuevo dataframe all.R podemos generar un gráfico donde se pueda observar el valor de \(\lambda\) que minimiza el \(\hat{R}(\lambda)\) el \(\hat{CV}(\lambda)\) y el \(\hat{GCV}(\lambda)\).
plot1 <- ggplot()+
geom_point(data = all.R, aes(x = LAMBDA, y = UBRE)) +
theme_bw() +
labs(x = expression(lambda), y = expression(hat(R)(lambda)))
plot2 <- ggplot()+
geom_point(data = all.R, aes(x = LAMBDA, y = CV)) +
theme_bw() +
labs(x = expression(lambda), y = expression(hat(CV)(lambda)))
plot3 <- ggplot()+
geom_point(data = all.R, aes(x = LAMBDA, y = GCV)) +
theme_bw() +
labs(x = expression(lambda), y = expression(hat(GCV)(lambda)))
grid.arrange(plot1, plot2, plot3, ncol=2)
tenemos entonces varios criterios y cada uno con una elección de \(\lambda\).
## Mejor lambda
## UBRE 9
## CV 9
## GCV 9## ubre cv GCV lambda
## 1 6.169403 60.45775 56.82861 9Observemos también el ajuste con \(\lambda=9\) con los datos reales (puntos) y los datos ajustados por el modelo (línea) de la variable \(NO_2\) vs \(hora\_normada\) para el día \(10\).
lambda <- 9
salida <- 2 # 1: Riesgo, 2: fitted values
fitted <- lambda.select(x, lambda, salida)
x <- bind_cols(x, fitted)ggplot()+
geom_point(data = x, aes(x = hora_normada, y = NO2)) +
geom_line(data = x, aes(x =hora_normada, y = fitted)) +
labs(y = expression(NO[2])) +
labs(title = "Estimación con base de cosenos") +
labs(subtitle = expression(lambda==9))
Tenemos entonces la estimación del comportamiento de \(NO_2\) para el dia \(10\) usando series de Fourier con base de cosenos y con un \(\lambda=9\), que seleccionamos por medio de los métodos \(UBRE\), \(CV\) y \(GCV\), podríamos decir que \(\hat{\mu}_9\) es una buena aproximación a \(\mu\).
El ajuste de este modelo se llevo acabo basándose en material expuesto en (Olaya, 2012).
Bibliografía
Olaya, J. (2012). Métodos de Regresión No Paramétrica. Universidad del Valle.
R Core Team. (2013). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. http://www.r-project.org/