ISLR实验：交叉验证法和自助法

February 10, 2021 (最后修改: November 08, 2024)

本文源自《统计学习导论：基于R语言应用》(ISLR) 中《5.3 实验：交叉验证法和自助法》章节

library(ISLR)

验证集方法

validation set approach

随机将观测集分成两部分：一个训练集 (training set) 和一个验证集 (validation set)，或叫做保留集 (hold-out set)。在训练集上拟合模型，用验证集计算错误率 (通常使用均方误差)。

head(Auto)

  mpg cylinders displacement horsepower weight acceleration year origin
1  18         8          307        130   3504         12.0   70      1
2  15         8          350        165   3693         11.5   70      1
3  18         8          318        150   3436         11.0   70      1
4  16         8          304        150   3433         12.0   70      1
5  17         8          302        140   3449         10.5   70      1
6  15         8          429        198   4341         10.0   70      1
                       name
1 chevrolet chevelle malibu
2         buick skylark 320
3        plymouth satellite
4             amc rebel sst
5               ford torino
6          ford galaxie 500

dim(Auto)

[1] 392   9

随机抽取训练集，使用 sample() 函数生成样本抽样序号子集。为保证结果一致，使用 set.seed() 函数设置固定的随机数种子。

set.seed(1)
train <- sample(nrow(Auto), nrow(Auto)/2)

拟合线性模型

lm_fit <- lm(
  mpg ~ horsepower,
  data=Auto,
  subset=train
)

计算验证集的均方误差

mean((Auto$mpg - predict(lm_fit, Auto))[-train]^2)

[1] 23.26601

计算二次和三次多项式回归的均方误差

lm_fit_p2 <- lm(
  mpg ~ poly(horsepower, 2),
  data=Auto,
  subset=train
)
mean((Auto$mpg - predict(lm_fit_p2, Auto))[-train]^2)

[1] 18.71646

lm_fit_p3 <- lm(
  mpg ~ poly(horsepower,3),
  data=Auto,
  subset=train
)
mean((Auto$mpg - predict(lm_fit_p3, Auto))[-train]^2)

[1] 18.79401

使用不同的训练集，在验证集上会得到不同的结果

set.seed(2)
train <- sample(nrow(Auto), nrow(Auto)/2)

lm_fit <- lm(
  mpg ~ horsepower,
  data=Auto,
  subset=train
)
mean((Auto$mpg - predict(lm_fit, Auto))[-train]^2)

[1] 25.72651

lm_fit_p2 <- lm(
  mpg ~ poly(horsepower,2),
  data=Auto,
  subset=train
)
mean((Auto$mpg - predict(lm_fit_p2, Auto))[-train]^2)

[1] 20.43036

lm_fit_p3 <- lm(
  mpg ~ poly(horsepower,3),
  data=Auto,
  subset=train
)
mean((Auto$mpg - predict(lm_fit_p3, Auto))[-train]^2)

[1] 20.38533

留一交叉验证法

leave-one-out cross-validation, LOOCV

将一个单独的观测作为验证集，剩下的观测作为训练集。多次重复该步骤，计算测试均方误差的 LOOCV 估计：

$$ CV_{(n)} = \frac {1} {n} \sum_{i=1}^n MSE_i $$

不设置 family 参数的 glm() 函数与 lm() 函数一样执行线性回归

glm_fit <- glm(
  mpg ~ horsepower,
  data=Auto
)
coefficients(glm_fit)

(Intercept)  horsepower 
 39.9358610  -0.1578447

lm_fit <- lm(
  mpg ~ horsepower,
  data=Auto
)
coefficients(lm_fit)

(Intercept)  horsepower 
 39.9358610  -0.1578447

glm() 与 boot 包的 cv.glm() 函数可以计算 LOOCV 估计

library(boot)

cv.glm() 函数返回结果中 delta 向量的两个数字为交叉验证的结果。本例中两个数字相等

glm_fit <- glm(
  mpg ~ horsepower,
  data=Auto
)
cv.err <- cv.glm(
  Auto,
  glm_fit
)
cv.err$delta

[1] 24.23151 24.23114

用递增的多项式次数重复上述步骤

cv.error <- rep(0, 5)
for (i in 1:5) {
  glm_fit <- glm(
    mpg ~ poly(horsepower, i),
    data=Auto
  )
  cv.error[i] <- cv.glm(
    Auto,
    glm_fit
  )$delta[1]
}
cv.error

[1] 24.23151 19.24821 19.33498 19.42443 19.03321

可以看到，超过二次的多项式拟合效果没有显著的提升

k 折交叉验证法

k-fold CV

将观测集随机地分为 k 个大小基本一致的组，或者说折 (fold)。将某一个折作为验证集，其余折作为训练集。重复该步骤 k 次，得到 k 折 CV 估计：

$$ CV_{(k)} = \frac {1} {k} \sum_{i=1}^k MSE_{i} $$

cv.glm() 函数也可以实现 k 折交叉验证法。通常 k 取 5 或 10。

set.seed(17)
cv_error_10 <- NULL
for (i in 1:10) {
  glm_fit <- glm(
    mpg ~ poly(horsepower, i),
    data=Auto
  )
  cv_error_10 <- rbind(
    cv_error_10,
    cv.glm(
      Auto,
      glm_fit,
      K=10
    )$delta
  )
}
cv_error_10

          [,1]     [,2]
 [1,] 24.27207 24.25459
 [2,] 19.26909 19.25393
 [3,] 19.34805 19.32647
 [4,] 19.29496 19.27197
 [5,] 19.03198 18.99974
 [6,] 18.89781 18.86178
 [7,] 19.12061 19.06318
 [8,] 19.14666 19.08645
 [9,] 18.87013 18.82284
[10,] 20.95520 20.74649

第一列是标准 k 折 CV 估计，第二列是偏差校正后的结果。本例中两者区别不大。

自助法

bootstrap

估计一个感兴趣的统计量的精度

金融资产数据集，X 和 Y 是两种金融资产的收益

head(Portfolio)

           X          Y
1 -0.8952509 -0.2349235
2 -1.5624543 -0.8851760
3 -0.4170899  0.2718880
4  1.0443557 -0.7341975
5 -0.3155684  0.8419834
6 -1.7371238 -2.0371910

dim(Portfolio)

[1] 100   2

计划将 alpha 比例投资 X，1 - alpha 比例投资到 Y。估计 alpha，使总风险或者说是总方差最小

alpha_fn <- function(data, index) {
  X <- data$X[index]
  Y <- data$Y[index]
  return(
    (var(Y) - cov(X, Y)) / (var(X) + var(Y) - 2*cov(X, Y))
  )
}

用全部数据估计

alpha_fn(Portfolio, 1:100)

[1] 0.5758321

sample() 有放回地取出 100 个观测进行估计 (replace=TRUE)

set.seed(1)
alpha_fn(Portfolio, sample(100, 100, replace=TRUE))

[1] 0.7368375

自助法，重复 1000 次上述步骤

boot_result <- boot(
  Portfolio,
  alpha_fn,
  R=1000
)
boot_result

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Portfolio, statistic = alpha_fn, R = 1000)


Bootstrap Statistics :
     original       bias    std. error
t1* 0.5758321 -0.001695873  0.09366347

估计 alpha 为 0.5758，标准误为 0.0937

使用 boot.ci() 函数求 alpha 的 95% 置信区间

boot.ci(boot_result, type="norm")

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 1000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot_result, type = "norm")

Intervals : 
Level      Normal        
95%   ( 0.3940,  0.7611 )  
Calculations and Intervals on Original Scale

估计线性回归模型的精度

衡量线性回归模型 beta_0 和 beta_1 估计的波动性

boot_fn <- function(data, index) {
  return (coef(lm(
    mpg ~ horsepower,
    data=Auto,
    subset=index
  )))
}

全部数据

boot_fn(Auto, 1:nrow(Auto))

(Intercept)  horsepower 
 39.9358610  -0.1578447

随机有放回的两个例子

set.seed(1)
boot_fn(
  Auto, 
  sample(
    nrow(Auto), 
    nrow(Auto), 
    replace=TRUE
  )
)
boot_fn(
  Auto, 
  sample(
    nrow(Auto), 
    nrow(Auto), 
    replace=TRUE
  )
)

(Intercept)  horsepower 
 40.3404517  -0.1634868 
(Intercept)  horsepower 
 40.1186906  -0.1577063

计算 1000 个截距和斜率

boot(
  Auto,
  boot_fn,
  R=1000
)

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Auto, statistic = boot_fn, R = 1000)


Bootstrap Statistics :
      original        bias    std. error
t1* 39.9358610  0.0544513229 0.841289790
t2* -0.1578447 -0.0006170901 0.007343073

用标准公式计算线性模型中截距和回归系数的标准误差

summary(
  lm(
    mpg ~ horsepower,
    data=Auto
  )
)$coef

              Estimate  Std. Error   t value      Pr(>|t|)
(Intercept) 39.9358610 0.717498656  55.65984 1.220362e-187
horsepower  -0.1578447 0.006445501 -24.48914  7.031989e-81

计算结果与自助法不同。实际上标准公式需要数据满足某些假设，而自助法则没有这些限制。数据有非线性关系，所以自助法得到的结果更接近真实值。

拟合二次模型的结果

boot_fn <- function(data, index) {
  return (coefficients(
    lm(
      mpg ~ horsepower + I(horsepower^2),
      data=data,
      subset=index
    )
  ))
}
set.seed(1)
boot(Auto, boot_fn, R=1000)

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Auto, statistic = boot_fn, R = 1000)


Bootstrap Statistics :
        original        bias     std. error
t1* 56.900099702  3.511640e-02 2.0300222526
t2* -0.466189630 -7.080834e-04 0.0324241984
t3*  0.001230536  2.840324e-06 0.0001172164

标准公式计算

summary(
  lm(
    mpg ~ horsepower + I(horsepower^2),
    data=Auto
  )
)$coef

                    Estimate   Std. Error   t value      Pr(>|t|)
(Intercept)     56.900099702 1.8004268063  31.60367 1.740911e-109
horsepower      -0.466189630 0.0311246171 -14.97816  2.289429e-40
I(horsepower^2)  0.001230536 0.0001220759  10.08009  2.196340e-21

二次模型拟合效果比线性模型更好，所以自助法估计与标准估计更接近

参考

https://github.com/perillaroc/islr-study

ISLR实验系列文章

线性回归

分类

重抽样方法

交叉验证法和自助法

线性模型选择与正则化