LOO-CV and GCV

本文展示了怎样使用LOO-CV与GCV（舍一交叉验证与广义交叉验证）交叉验证方法去寻找LOESS（局部加权多项式回归）中的最优span并展示出最优span下真实曲线与回归曲线的差异。

LOO-CV and GCV

LOO-CV prediction error can be computed based on the original model (the one that uses all the n samples):

\[\frac{1}{n}\sum_{i=1}^n (y_i-\hat{y_i}^{[-i]})^2=\frac{1}{n}\sum_{i=1}^n(\frac{y_i-\hat{y_i}}{1-S_{ii}})^2\] While \(S_{ii}\) is the (i, i)-th entry of the smoothing matrix S

GCV prediction error is very similar to LOO-CV

\[\frac{1}{n}\sum_{i=1}^n(\frac{y_i-\hat{y_i}}{1-trace(S)})^2\]

library(tidyverse)

## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --

## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.4     v dplyr   1.0.7
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   2.0.1     v forcats 0.5.1

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

lo.lev <- function(x1, sp){
  # x1: n-by-1 feature vector
  # sp: a numerical value for "span"
  
  n = length(x1);
  lev = rep(0, n)
  
  
  for(i in 1:n){
       y = rep(0, n); y[i]=1;
       si = loess(y ~ x1, sp=sp)$fitted
       si=si[i]
       lev[i]= si;
       }
  
  return(lev)
}

onestep_CV <- function(x1, y1, sp){
  
  CVmode=loess(y1 ~ x1, sp=sp)
  lolev=lo.lev(x1,sp)
  m=mean(lolev)
  cv=mean(((y1-CVmode$fitted)/(1-lolev))^2)
  gcv=mean(((y1-CVmode$fitted)/(1-m))^2)
  
  return(list(cv = cv, gcv = gcv))
}

myCV <- function(x1, y1, span){
  # x1: feature vector of length n
  # y1: response vector of length n
  # span: a sequence of values for "span"
  
  m = length(span)
  cv = rep(0, m)
  gcv = rep(0, m)
  
  for(i in 1:m){
    tmp = onestep_CV(x1, y1, span[i])
    cv[i] = tmp$cv
    gcv[i] = tmp$gcv
  }
  return(list(cv = cv, gcv = gcv))
}

The three functions above are the key function to obtain the estimation error for different spans range from 0.2 to 0.9.

lo.lev function is to calculate the (i, i)-th entry of the smoothing matrix of LOESS, onestep_CV is to find the estimation error given by LOO-CV and GCV with a fixed span. The main function myCV is to give a loop to find the error under different spans.

And then generating data

mydata=read_csv("https://raw.githubusercontent.com/Zhenz2020/Zhenz2020.github.io/main/Data_LOOCV.csv")

## Rows: 30 Columns: 2

## -- Column specification --------------------------------------------------------
## Delimiter: ","
## dbl (2): x, y

## 
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.

span_value=seq(from=0.2,by=0.05,length=15)

cv.out = myCV(mydata$x, mydata$y, span_value)
myout = data.frame(CV = cv.out$cv, 
                   GCV = cv.out$gcv, 
                   span = span_value)
myout

##           CV      GCV span
## 1  12.415911 2.110162 0.20
## 2   2.241473 1.489206 0.25
## 3   1.502980 1.190110 0.30
## 4   1.259175 1.174423 0.35
## 5   1.190380 1.102540 0.40
## 6   1.156812 1.062503 0.45
## 7   1.124306 1.040422 0.50
## 8   1.179664 1.118841 0.55
## 9   1.179464 1.119269 0.60
## 10  1.250914 1.180585 0.65
## 11  1.553562 1.519091 0.70
## 12  1.636175 1.627429 0.75
## 13  1.764534 1.744549 0.80
## 14  1.976094 1.925696 0.85
## 15  2.035108 1.979820 0.90

myout$span[myout$GCV == min(myout$GCV)]

## [1] 0.5

myout$span[myout$CV == min(myout$CV)]

## [1] 0.5

As we can see above, the best spans for LOO-CV and GCV are both 0.5

And then, we can plot the the data and the estimation curve by span=0.5

spangcv.min = 0.5
plot(mydata$x, mydata$y, xlab="", ylab="", col="gray");
fx = 1:50/50;
fy = sin(12*(fx+0.2))/(fx+0.2)
lines(fx, fy, col=8, lwd=2);
f = loess(y ~ x, mydata, span = spangcv.min)
lines(fx, predict(f, data.frame(x = fx), surface = "direct"), 
      lty=2, lwd=2, col="blue")

In the figure above gray circle is the data, true curve is the gray one and the LOOES curve is the blue one.