모델의 성능을 평가할 때에는 항상 검증셋(test dataset)을 이용해야 합니다. training set에서는 항상 더 유연한 모델이 더 정확한 예측 결과를 보여줍니다. 그러나, 훈련셋에서 더 정확한 결과가 검정셋에서의 정확성을 보장해주지는 못합니다. 오히려, 유연한 모델일 수록 검정셋의 noise까지 학습하게 됩니다.즉, overfit 하게 되어 variance 가 커지게 됩니다. 우리는 bias-variance trade-off 를 고려하여 적절한 수준의 flexibility 를 갖는 모델을 선택(model selection)해야 합니다.
검증셋을 이용하여 모델의 성능을 평가하는 방법에는 아래의 세가지가 있습니다.
Auto 데이터셋에서 연비(mpg)를 예측하는 회귀 모델을 만들고, 모델의 성능을 평가해봅시다.
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
Auto = read.csv("Auto.csv")
str(Auto)
## 'data.frame': 397 obs. of 9 variables:
## $ mpg : num 18 15 18 16 17 15 14 14 14 15 ...
## $ cylinders : int 8 8 8 8 8 8 8 8 8 8 ...
## $ displacement: num 307 350 318 304 302 429 454 440 455 390 ...
## $ horsepower : Factor w/ 94 levels "?","100","102",..: 17 35 29 29 24 42 47 46 48 40 ...
## $ weight : int 3504 3693 3436 3433 3449 4341 4354 4312 4425 3850 ...
## $ acceleration: num 12 11.5 11 12 10.5 10 9 8.5 10 8.5 ...
## $ year : int 70 70 70 70 70 70 70 70 70 70 ...
## $ origin : int 1 1 1 1 1 1 1 1 1 1 ...
## $ name : Factor w/ 304 levels "amc ambassador brougham",..: 49 36 231 14 161 141 54 223 241 2 ...
Auto$horsepower = as.numeric(as.character(Auto$horsepower))
## Warning: NAs introduced by coercion
Auto = Auto[complete.cases(Auto),]
summary(Auto)
## mpg cylinders displacement horsepower
## Min. : 9.00 Min. :3.000 Min. : 68.0 Min. : 46.0
## 1st Qu.:17.00 1st Qu.:4.000 1st Qu.:105.0 1st Qu.: 75.0
## Median :22.75 Median :4.000 Median :151.0 Median : 93.5
## Mean :23.45 Mean :5.472 Mean :194.4 Mean :104.5
## 3rd Qu.:29.00 3rd Qu.:8.000 3rd Qu.:275.8 3rd Qu.:126.0
## Max. :46.60 Max. :8.000 Max. :455.0 Max. :230.0
##
## weight acceleration year origin
## Min. :1613 Min. : 8.00 Min. :70.00 Min. :1.000
## 1st Qu.:2225 1st Qu.:13.78 1st Qu.:73.00 1st Qu.:1.000
## Median :2804 Median :15.50 Median :76.00 Median :1.000
## Mean :2978 Mean :15.54 Mean :75.98 Mean :1.577
## 3rd Qu.:3615 3rd Qu.:17.02 3rd Qu.:79.00 3rd Qu.:2.000
## Max. :5140 Max. :24.80 Max. :82.00 Max. :3.000
##
## name
## amc matador : 5
## ford pinto : 5
## toyota corolla : 5
## amc gremlin : 4
## amc hornet : 4
## chevrolet chevette: 4
## (Other) :365
set.seed(1)
train_index = sample(nrow(Auto), round(nrow(Auto)*2/3))
Auto = select(Auto, -name)
train_data = Auto[train_index, ]
test_data = Auto[-train_index, ]
lm.fit = lm(mpg ~ ., data = train_data)
pred = predict(lm.fit, test_data)
obs = test_data$mpg
mse = mean((pred - obs)^2)
mse
## [1] 10.56468
validation set approach 를 10번 반복해서 평균제곱오차의 추정치와 분산을 구해보자.
mse_vset = c()
for (i in 1:10){
train_index = sample(nrow(Auto), round(nrow(Auto)*2/3), replace = F)
train_data = Auto[train_index, ]
test_data = Auto[-train_index, ]
lm.fit = lm(mpg ~ ., data = train_data)
pred = predict(lm.fit, test_data)
obs = test_data$mpg
mse_vset[i] = mean((pred - obs)^2)
}
mean(mse_vset)
## [1] 11.35968
var(mse_vset)
## [1] 3.006118
LOOCV 를 이용하여 평균제곱오차의 추정치와 분산을 구해보자.
mse_loocv = c()
for (i in 1:nrow(Auto)){
lm.fit = lm(mpg ~ ., data = Auto[-i,])
pred = predict(lm.fit, Auto[i,])
obs = Auto$mpg[i]
mse_loocv[i] = (pred - obs)^2
}
mean(mse_loocv)
## [1] 11.37113
var(mse_loocv)
## [1] 447.8021
K-fold CV
library(caret)
## Loading required package: lattice
## Loading required package: ggplot2
set.seed(1)
folds = createFolds(1:nrow(Auto), k=10)
mse_kcv = c()
for (i in 1:10){
lm.fit = lm(mpg ~ ., data = Auto[-folds[[i]],])
pred = predict(lm.fit, Auto[folds[[i]],])
obs = Auto[folds[[i]],]$mpg
mse_kcv[i] = mean((pred - obs)^2)
}
mean(mse_kcv)
## [1] 11.34605
var(mse_kcv)
## [1] 8.582353
위 세가지 검증 기법을 비교해봅시다.
Validation set approach는 training set의 크기가 적어 bias가 큰 경향이 있고, training vs. test dataset 분할의 임의성으로 인해 test error 에 대한 추정치의 변동성이 보통 크게 나타납니다. 즉, bias와 variance가 모두 큽니다.
반면에, LOOCV는 training dataset의 크기가 커서 bias를 줄일 수 있고, k-fold CV 에서의 bias 는 validation set approach와 LOOCV의 사이에 해당합니다.
variance 는 LOOCV 와 k-fold CV 중 어느 것이 더 클까요? 보통 예상과 달리 k-fold CV 이 LOOCV 보다 variance 가 더 작은 경향이 있습니다. LOOCV 는 training dataset 이 서로 매우 유사하므로 모델간의 상관 관계가 매우 높습니다. 모델 간의 상관 관계가 높을 수록 variance 가 더 크게 나타납니다.
추정치의 정확도, 즉 변동성 평가
표본평균에서 모평균을 추정하는 문제를 기억해봅시다.
다음과 같은 모집단이 있다고 합시다.
population = rnorm(1000, 0, 10) # number=1000, mean=0, sd=10
mean(population)
## [1] -0.4106525
sd(population)
## [1] 10.51136
모집단에서 표본을 추출합니다.
sampl = sample(population, 100, replace = F)
이제, 모집단의 평균과 표준편차를 모른다고 가정하고, 표본 데이터를 이용해서 모평균과 모표준편차를 추정해봅시다.
먼저, 모평균을 추청해봅시다.
중심극한정리에 따라서 모집단의 평균은 표본평균의 평균과 같습니다.
\[\mu = mean(\bar{X})\] \(\mu: population\ mean\)
\(\bar{X}: sample\ mean\)
즉, 표본평균의 평균이 모평균의 추정치가 되는데, 이 추정치는 얼마나 정확할까요? 우리는 먼저 모평균을 모른다고 가정하였으므로 추정치가 얼마나 정확한지 알수 없습니다. 따라서, 추정치가 얼마나 정확한가에 대한 질문보다는 추정치의 변동성이 얼마나 큰가라는 질문이 더 적절한 질문입니다. 추정치의 변동성을 나타내는 지표가 바로 표준 오차입니다.
standard error (of mean) = standard deviation of estimated population mean \[SE=\frac{\sigma}{\sqrt{n}}\]
\(\sigma = standard\ deviation\ of\ population\) \(n = sample\ size\)
참고로, 신뢰구간은 표준오차와 비슷한 개념인데(추정치의 변동성) 모집단의 모수값이 포함될 가능성이 있는 범위를 나타냅니다.
95% confidence interval of population mean \[\bar{X}-1.96\times\frac{\sigma}{\sqrt{n}} < \mu < \bar{X}+1.96\times\frac{\sigma}{\sqrt{n}}\]
다시, 표준오차로 돌아가서 우리는 모집단의 표준편차 \(\sigma\)를 모른다고 가정하였으므로, 표본의 표준편차 \(\hat{\sigma}\)를 대신 사용합니다.
\[\hat{\sigma} = s\] \(s = standard\ deviation\ of\ sample\)
표본의 표준편차는 다음과 같이 구할 수 있습니다.
\[s^2 = \frac{\sum(X_i-\bar{X})^2}{n-1}\] note: n-1 instead of n 왜, n 대신 n-1을 사용하는 걸까요?
위 식과 같이 n-1을 사용해서 구한 표준편차는 정확히 말하면 표본 표준편차는 아닙니다. 대신 우리는 이를 모집단의 표준편차에 대한 불편 추정량이라고 합니다. 모집단의 표준편차는 항상 표본의 표준편차보다 클 것입니다. 1을 빼주는 것은 우리는 표본의 평균을 알고 있고, 표본의 평균이 주어지면 잔차의 합은 항상 0이기 때문입니다. 즉, 모집단의 분산에 대한 불편 추정량으로서 표본의 분산에 대한 자유도는 n이 아니라 n-1 이 됩니다.
이제 위에서 설명한 식을 이용해서 다음과 같이 표본에서 모평균을 추정할 수 있습니다.
x_bar = mean(sampl)
n = length(sampl)
s = sqrt((sd(sampl)^2)*n/(n-1))
se = s/sqrt(n)
CI = c(x_bar - 1.96*s/sqrt(n), x_bar + 1.96*s/sqrt(n))
se; CI
## [1] 1.146679
## [1] -2.271601 2.223380
다음으로 bootstrapping 기법을 이용해 모평균을 추정해봅시다. bootstrapping 을 위해 반복적으로 적용한 function 을 먼저 정의해주어야 합니다. 여기서는 표본의 평균을 추정하는 함수가 되겠습니다.
fn = function(z, index){
mean(z[index])
} # function need two arguments, the second one should be index
이제 표본에서 하위 표본을 추출하고(복원을 허락, resample), 해당 하위 표본에 평균을 구하는 함수를 적용합니다.
fn(sampl, sample(length(sampl), length(sampl), replace = T))
## [1] 0.09466536
위 과정을 n번 반복합니다.
# boot
library(boot)
##
## Attaching package: 'boot'
## The following object is masked from 'package:lattice':
##
## melanoma
res = boot(sampl, fn, R=1000)
res$t0
## [1] -0.02411041
# original = the mean of the original sample
mean(res$t)
## [1] -0.00940608
mean(res$t) - res$t0
## [1] 0.01470433
# bias = the mean of resamples (the esimated mean of the population) - the mean of original sample
sd(res$t)
## [1] 1.169339
# std. error = the standard deviation of the mean of resamples
boot.ci(res)
## Warning in boot.ci(res): bootstrap variances needed for studentized
## intervals
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = res)
##
## Intervals :
## Level Normal Basic
## 95% (-2.3307, 2.2530 ) (-2.4767, 2.2449 )
##
## Level Percentile BCa
## 95% (-2.2931, 2.4285 ) (-2.2130, 2.5167 )
## Calculations and Intervals on Original Scale
bootstrapping 방법으로 추정한 모평균과 신뢰구간을 위에서 모수적으로 추정한 모평균, 신뢰구간과 비교해봅시다.