# Cross-validation and bootstep
# Source: Course 'Math modeling' practical work, State University of Management, Moscow
# link: https://sites.google.com/a/kiber-guu.ru/r-practice/home
library('ISLR') # набор данных Auto
library('GGally') # матричные графики
library('boot') # расчёт ошибки с кросс-валидацией
my.seed <- 1
# Get data and do primary visual inspection
head(Auto)
str(Auto)
ggpairs(Auto[, -9])
# только mpg ~ horsepower
plot(Auto$horsepower, Auto$mpg,
xlab = 'horsepower', ylab = 'mpg', pch = 21,
col = rgb(0, 0, 1, alpha = 0.4), bg = rgb(0, 0, 1, alpha = 0.4))
# общее число наблюдений
n <- nrow(Auto)
# доля обучающей выборки
train.percent <- 0.5
# выбрать наблюдения в обучающую выборку
set.seed(my.seed)
inTrain <- sample(n, n * train.percent)
plot(Auto$horsepower[inTrain], Auto$mpg[inTrain],
xlab = 'horsepower', ylab = 'mpg', pch = 21,
col = rgb(0, 0, 1, alpha = 0.4), bg = rgb(0, 0, 1, alpha = 0.4))
points(Auto$horsepower[-inTrain], Auto$mpg[-inTrain],
pch = 21, col = rgb(1, 0, 0, alpha = 0.4), bg = rgb(1, 0, 0, alpha = 0.4))
legend('topright',
pch = c(16, 16), col = c('blue', 'red'), legend = c('test', 'train'))
### Build different types of models
## 1. Linear model mpg = a+b*horsepower
# присоединить таблицу с данными: названия стоблцов будут доступны напрямую
attach(Auto)
# подгонка линейной модели на обучающей выборке
fit.lm.1 <- lm(mpg ~ horsepower, subset = inTrain)
# считаем MSE на тестовой выборке
mean((mpg[-inTrain] - predict(fit.lm.1, Auto[-inTrain, ]))^2)
# отсоединить таблицу с данными
detach(Auto)
#------------------------------------------------------------------------------
## 2. Square Model mpg = a + b*horsepower + c*horsepower^2
# присоединить таблицу с данными: названия стоблцов будут доступны напрямую
attach(Auto)
# подгонка линейной модели на обучающей выборке
fit.lm.2 <- lm(mpg ~ poly(horsepower, 2), subset = inTrain)
# считаем MSE на тестовой выборке
mean((mpg[-inTrain] - predict(fit.lm.2, Auto[-inTrain, ]))^2)
detach(Auto)
#------------------------------------------------------------------------------
## 3. Qubiq Model mpg = a + b*horsepower + c*horsepower^2 + d*horsepower^3
# присоединить таблицу с данными: названия стоблцов будут доступны напрямую
attach(Auto)
# подгонка линейной модели на обучающей выборке
fit.lm.3 <- lm(mpg ~ poly(horsepower, 3), subset = inTrain)
# считаем MSE на тестовой выборке
mean((mpg[-inTrain] - predict(fit.lm.3, Auto[-inTrain, ]))^2)
detach(Auto)
#------------------------------------------------------------------------------
## LOOCV - Cross validation
# подгонка линейной модели на обучающей выборке
fit.glm <- glm(mpg ~ horsepower, data = Auto)
# считаем LOOCV-ошибку
cv.err <- cv.glm(Auto, fit.glm)
# результат: первое число -- по формуле LOOCV-ошибки,
# второе -- с поправкой на смещение
cv.err$delta[1]
# вектор с LOOCV-ошибками
cv.err.loocv <- rep(0, 5)
names(cv.err.loocv) <- 1:5
# цикл по степеням полиномов
for (i in 1:5){
fit.glm <- glm(mpg ~ poly(horsepower, i), data = Auto)
cv.err.loocv[i] <- cv.glm(Auto, fit.glm)$delta[1]
}
# результат
cv.err.loocv
#------------------------------------------------------------------------------
## k-something cross-validation k-кратная перекрёстная проверка
# оценим точность полиномиальных моделей, меняя степень
# вектор с ошибками по 10-кратной кросс-валидации
cv.err.k.fold <- rep(0, 5)
names(cv.err.k.fold) <- 1:5
# цикл по степеням полиномов
for (i in 1:5){
fit.glm <- glm(mpg ~ poly(horsepower, i), data = Auto)
cv.err.k.fold[i] <- cv.glm(Auto, fit.glm,
K = 10)$delta[1]
}
# result
cv.err.k.fold
# Compare with MSE
err.test
##### BOOTSTREP
head(Portfolio)
str(Portfolio)
# функция для вычисления искомого параметра
alpha.fn <- function(data, index){
X = data$X[index]
Y = data$Y[index]
(var(Y) - cov(X, Y)) / (var(X) + var(Y) - 2*cov(X, Y))
}
# рассчитать alpha по всем 100 наблюдениям
alpha.fn(Portfolio, 1:100)
# создать бутстреп-выборку и повторно вычислить alpha
set.seed(my.seed)
alpha.fn(Portfolio, sample(100, 100, replace = T))
# теперь -- многократное повторение предыдущей операции
boot(Portfolio, alpha.fn, R = 1000)
#Бутстреп повторяет расчёт параметра много раз, делая повторные выборки из наших 100 наблюдений. В итоге этим методом можно вычислить стандартную ошибку параметра, не опираясь на допущения о законе распределении параметра. В нашем случае α=0.576 со стандартной ошибкой sα^=0.089.
#Tочность оценки параметра регрессии
# При построении модели регрессии проблемы в остатках приводят к неверной оценке ошибок параметров. Обойти эту проблему можно, применив для расчёта этих ошибок бутстреп.
# Оценивание точности линейной регрессионной модели ----------------------------
# оценить стандартные ошибки параметров модели
# mpg = beta_0 + beta_1 * horsepower с помощью бутстрепа,
# сравнить с оценками ошибок по МНК
# функция для расчёта коэффициентов ПЛР по выборке из данных
boot.fn <- function(data, index){
coef(lm(mpg ~ horsepower, data = data, subset = index))
}
boot.fn(Auto, 1:n)
# пример применения функции к бутстреп-выборке
set.seed(my.seed)
boot.fn(Auto, sample(n, n, replace = T))
# применяем функцию boot для вычисления стандартных ошибок параметров
# (1000 выборок с повторами)
boot(Auto, boot.fn, 1000)
# сравним с МНК
attach(Auto)
summary(lm(mpg ~ horsepower))$coef
detach(Auto)
# оценки отличаются из-за того, что МНК -- параметрический метод с допущениями
# вычислим оценки параметров квадратичной модели регрессии
boot.fn.2 <- function(data, index){
coef(lm(mpg ~ poly(horsepower, 2), data = data, subset = index))
}
# применим функцию к 1000 бутсреп-выборкам
set.seed(my.seed)
boot(Auto, boot.fn, 1000)