Практика 1
Математическое моделирование
Оценка точности модели с непрерывной зависимой переменной
В практических примерах ниже показано:
- как делить данные на выборки (обучающую и тестовую);
- как считать MSE: среднеквадратическую ошибку модели;
- как меняются MSE на тестовой и обучающей выборках с изменением гибкости (числа степеней свободы) модели.
Модели: сглаживающие сплайны.
Данные: сгенерированные.
Рассмотрим пример из лекции: как меняется поведение ошибок на тестовой и обучающей выборках при различном числе степеней свободы, если функция зависимости отклика \(Y\) от единственного признака \(X\) известна. Сгенерируем \(X\) и \(Y\):
- \(X \sim U(5, 105)\)
- \(Y = f(X) + \epsilon\), где \(f(X) = 4 - 0.02X + 0.0055X^2 - 4.9 \cdot 10^{-5} \cdot X^3\); \(\epsilon \sim N(0, 1)\).
# Генерируем данные ###########################################################
my.seed <- 1486372882 # ядро
n.all <- 60 # наблюдений всего
train.percent <- 0.85 # доля обучающей выборки
res.sd <- 1 # стандартное отклонение случайного шума
x.min <- 5 # границы изменения X: нижняя
x.max <- 105 # и верхняя
# фактические значения x
set.seed(my.seed)
x <- runif(x.min, x.max, n = n.all)
# случайный шум
set.seed(my.seed)
res <- rnorm(mean = 0, sd = res.sd, n = n.all)
# отбираем наблюдения в обучающую выборку
set.seed(my.seed)
inTrain <- sample(seq_along(x), size = train.percent*n.all)
# истинная функция взаимосвязи
y.func <- function(x) {4 - 2e-02*x + 5.5e-03*x^2 - 4.9e-05*x^3}
# для графика истинной взаимосвязи
x.line <- seq(x.min, x.max, length = n.all)
y.line <- y.func(x.line)
# фактические значения y (с шумом)
y <- y.func(x) + res
# Создаём векторы с данными для построения графиков ############################
# наблюдения на обучающей выборке
x.train <- x[inTrain]
y.train <- y[inTrain]
# наблюдения на тестовой выборке
x.test <- x[-inTrain]
y.test <- y[-inTrain]Изобразим исходные данные на графике.
# График 1: Исходные данные на график #########################################
# убираем широкие поля рисунка
par(mar = c(4, 4, 1, 1))
# наименьшие/наибольшие значения по осям
x.lim <- c(x.min, x.max)
y.lim <- c(min(y), max(y))
# наблюдения с шумом (обучающая выборка)
plot(x.train, y.train,
col = grey(0.2), bg = grey(0.2), pch = 21,
xlab = 'X', ylab = 'Y',
xlim = x.lim, ylim = y.lim,
cex = 1.2, cex.lab = 1.2, cex.axis = 1.2)
# заголовок
mtext('Исходные данные и истинная функция связи', side = 3)
# наблюдения тестовой выборки
points(x.test, y.test, col = 'red', bg = 'red', pch = 21)
# истинная функция
lines(x.line, y.line, lwd = 2, lty = 2)
# легенда
legend('topleft', legend = c('обучение', 'тест', 'f(X)'),
pch = c(16, 16, NA),
col = c(grey(0.2), 'red', 'black'),
lty = c(0, 0, 2), lwd = c(1, 1, 2), cex = 1.2)
В качестве модели используем сплайны со степенями свободы от 2 (прямая) до 40 (количество узлов равно 2/3 наблюдений). Строим модели с различным количеством степеней свободы и в каждом случае считаем среднеквадратическую ошибку модели на обучающей и тестовой выборках.
# Строим модель №2 из лекции (df = 6) #########################################
# модель 2 (сплайн с df = 6)
mod <- smooth.spline(x = x.train, y = y.train, df = 6)
# модельные значения для расчёта ошибок
y.model.train <- predict(mod, data.frame(x = x.train))$y[, 1]
y.model.test <- predict(mod, data.frame(x = x.test))$y[, 1]
# считаем средний квадрат ошибки на обечающей и тестовой выборке
MSE <- c(sum((y.train - y.model.train)^2) / length(x.train),
sum((y.test - y.model.test)^2) / length(x.test))
names(MSE) <- c('train', 'test')
round(MSE, 2)## train test
## 0.95 1.21# Теперь строим модели с df от 2 до 40 ########################################
# максимальное число степеней свободы для модели сплайна
max.df <- 40
tbl <- data.frame(df = 2:max.df) # таблица для записи ошибок
tbl$MSE.train <- 0 # столбец: ошибки на обучающей выборке
tbl$MSE.test <- 0 # столбец: ошибки на тестовой выборке
# цикл по степеням свободы
for (i in 2:max.df) {
# строим модель
mod <- smooth.spline(x = x.train, y = y.train, df = i)
# модельные значения для расчёта ошибок
y.model.train <- predict(mod, data.frame(x = x.train))$y[, 1]
y.model.test <- predict(mod, data.frame(x = x.test))$y[, 1]
# считаем средний квадрат ошибки на обучающей и тестовой выборке
MSE <- c(sum((y.train - y.model.train)^2) / length(x.train),
sum((y.test - y.model.test)^2) / length(x.test))
# записываем ошибки в таблицу
tbl[tbl$df == i, c('MSE.train', 'MSE.test')] <- MSE
}
# первые строки таблицы
head(tbl)Изобразим на графике поведение ошибок при различном количестве степеней свободы.
# График 2: Зависимость MSE от гибкости модели ################################
plot(x = tbl$df, y = tbl$MSE.test,
type = 'l', col = 'red', lwd = 2,
xlab = 'Степени свободы сплайна', ylab = 'MSE',
ylim = c(min(tbl$MSE.train, tbl$MSE.test),
max(tbl$MSE.train, tbl$MSE.test)),
cex = 1.2, cex.lab = 1.2, cex.axis = 1.2)
# заголовок
mtext('Изменение MSE с ростом числа степеней свободы', side = 3)
points(x = tbl$df, y = tbl$MSE.test,
pch = 21, col = 'red', bg = 'red')
lines(x = tbl$df, y = tbl$MSE.train, col = grey(0.3), lwd = 2)
# неустранимая ошибка
abline(h = res.sd, lty = 2, col = grey(0.4), lwd = 2)
# легенда
legend('topleft', legend = c('обучающая', 'тестовая'),
pch = c(NA, 16),
col = c(grey(0.2), 'red'),
lty = c(1, 1), lwd = c(2, 2), cex = 1.2)
# степени свободы у наименьшей ошибки на тестовой выборке
min.MSE.test <- min(tbl$MSE.test)
df.min.MSE.test <- tbl[tbl$MSE.test == min.MSE.test, 'df']
# компромисс между точностью и простотой модели по графику
df.my.MSE.test <- 6
my.MSE.test <- tbl[tbl$df == df.my.MSE.test, 'MSE.test']
# ставим точку на графике
abline(v = df.my.MSE.test,
lty = 2, lwd = 2)
points(x = df.my.MSE.test, y = my.MSE.test,
pch = 15, col = 'blue')
mtext(df.my.MSE.test,
side = 1, line = -1, at = df.my.MSE.test, col = 'blue', cex = 1.2)
На этом графике:
- При движении слева направо MSE на обучающей выборке (серая кривая) сокращается, потому что с ростом числа степеней свободы расчёт число узлов, по которым строится сплайн. При этом модельная кривая подгоняется по всё возрастающему количеству точек и становится всё более гибкой. В результате индивидуальные расстояния от фактических наблюдений за \(Y\) до их модельных оценок сокращаются, что приводит к сокращению MSE.
- При движении слева направо MSE на тестовой выборке (красная кривая) сначала резко сокращается, затем растёт. Нам известна истинная форма связи \(Y\) с \(X\), она описывается кубической функцией. Число степеней свободы такой модели равно числу оцениваемых параметров, т.е. 4 (коэффициенты перед \(X\), \(X^2\), \(X^3\) и константа). Поэтому резкое падение ошибки на тестовой выборке при небольшом числе степеней свободы связано с тем, что модель приближается по гибкости к истинной функции связи. Затем MSE на тестовой выборке довольно долго остаётся стабильной, а затем начинает расти. Этот рост объясняется эффектом переобучения модели: она всё лучше описывает обучающую выборку, и при этом постепенно становится неприменимой ни к одному другому набору наблюдений.
Наименьшее значение MSE на тестовой выборке соответствует числу степеней свободы 7 и равно 1.21. Визуально по графику мы можем установить, что первое значение \(MSE_{ТЕСТ}\), близкое к стабильно низким, соответствует df = 6. Ошибка здесь равна 1, что ненамного отличается от минимума. Именно df = 6 было выбрано в качестве компромисса между точностью (минимальной MSE на тестовой выборке) и простотой модели (чем меньше степеней свободы, тем модель проще).
График с моделью, выбранной в качестве лучшей, показан на рисунке ниже.
# График 3: Лучшая модель (компромисс между гибкостью и точностью) ############
mod.MSE.test <- smooth.spline(x = x.train, y = y.train, df = df.my.MSE.test)
# для гладких графиков модели
x.model.plot <- seq(x.min, x.max, length = 250)
y.model.plot <- predict(mod.MSE.test, data.frame(x = x.model.plot))$y[, 1]
# убираем широкие поля рисунка
par(mar = c(4, 4, 1, 1))
# наименьшие/наибольшие значения по осям
x.lim <- c(x.min, x.max)
y.lim <- c(min(y), max(y))
# наблюдения с шумом (обучающая выборка)
plot(x.train, y.train,
col = grey(0.2), bg = grey(0.2), pch = 21,
xlab = 'X', ylab = 'Y',
xlim = x.lim, ylim = y.lim,
cex = 1.2, cex.lab = 1.2, cex.axis = 1.2)
# заголовок
mtext('Исходные данные и лучшая модель', side = 3)
# наблюдения тестовой выборки
points(x.test, y.test,
col = 'red', bg = 'red', pch = 21)
# истинная функция
lines(x.line, y.line,
lwd = 2, lty = 2)
# модель
lines(x.model.plot, y.model.plot,
lwd = 2, col = 'blue')
# легенда
legend('topleft', legend = c('обучение', 'тест', 'f(X)', 'модель'),
pch = c(16, 16, NA, NA),
col = c(grey(0.2), 'red', 'black', 'blue'),
lty = c(0, 0, 2, 1), lwd = c(1, 1, 2, 2), cex = 1.2)
Упражнение 1
Завести аккаунт на github.com.
Задача 1. На данных своего варианта повторить три графика из первой практики, выбрав число степеней свободы как компромисс между точностью (оценкой ошибки на тестовой выборке) и простотой модели (числом степеней свободы). Все рисунки сохранить в графические файлы в формате png.
Задача 2. Решить задачу 1, изменив характеристики данных (см. свой вариант). Почему при таком изменении данных MSE меняется именно так? Все рисунки сохранить в графические файлы в формате png.
- Выполненные задачи 1-2 разместить в одном отчёте в репозитории на github.com, выслать ссылку на него на почту преподавателя. В репозитории должны лежать:
- файл README.md с кратким описанием содержимого репозитория и с ответами на вопросы задач;
- скрипт генерации отчёта: файл .Rmd в кодировке UTF-8;
- графики, сгенерированные в задании.
В отчёте с решением должны присутствовать, кроме блоков кода R, вводный текст с постановкой задачи и анализ динамики MSE на тестовой и обучающей выборках в вашем варианте.
Варианты
Все условия, не упомянутые в таблице (величина выборки, закон распределения \(X\) и т.д.) брать из примеров первой практики.
Номер варианта – номер студента в списке. Студент под номером 21 берёт вариант 1, под номером 22 – 2, и т.д.
В качестве ядра генератора случайных чисел используйте номер своего варианта.
| Номер варианта | Функция для задачи 1 | Характеристики для задачи 2 |
| 1 | \(f(X) = 13 + 3.5 \cdot \sin{\Big ( {x - 30 \over 9 } \Big )}\) | \(n = 600\), \(n = 550\), \(n = 500\) |
| 2 | \(f(X) = 17 + 0.02 \cdot x -0.005 \cdot (x - 45)^2 + 0.00006 \cdot (x - 54)^3\) | \(train.percent = 0.9\), \(train.percent = 0.85\), \(train.percent = 0.8\) |
| 3 | \(f(X) = 19 - 0.05 \cdot x\) | \(res.sd = 2\), \(res.sd = 3\), \(res.sd = 4\) |
| 4 | \(f(X) = 11 + 3.5 \cdot \sin{\Big ( {x - 30 \over 9 } \Big )}\) | \(train.percent = 0.3\), \(train.percent = 0.25\), \(train.percent = 0.2\) |
| 5 | \(f(X) = 15 + 0.02 \cdot x -0.005 \cdot (x - 45)^2 + 0.00006 \cdot (x - 54)^3\) | \(res.sd = 0.5\), \(res.sd = 1\), \(res.sd = 1.5\) |
| 6 | \(f(X) = 21 - 0.05 \cdot x\) | \(n = 450\), \(n = 400\), \(n = 350\) |
| 7 | \(f(X) = 5 + 3.5 \cdot \sin{\Big ( {x - 30 \over 9 } \Big )}\) | \(res.sd = 5.5\), \(res.sd = 6\), \(res.sd = 6.5\) |
| 8 | \(f(X) = 15 + 0.02 \cdot x -0.005 \cdot (x - 45)^2 + 0.00006 \cdot (x - 54)^3\) | \(n = 300\), \(n = 250\), \(n = 200\) |
| 9 | \(f(X) = 20 - 0.05 \cdot x\) | \(train.percent = 0.5\), \(train.percent = 0.45\), \(train.percent = 0.4\) |
| 10 | \(f(X) = 6 + 3.5 \cdot \cos{\Big ( {x - 30 \over 12 } \Big )}\) | \(n = 150\), \(n = 100\), \(n = 50\) |
| 11 | \(f(X) = 25 + 0.02 \cdot x -0.003 \cdot (x - 45)^2 + 0.00006 \cdot (x - 54)^3\) | \(train.percent = 0.2\), \(train.percent = 0.15\), \(train.percent = 0.1\) |
| 12 | \(f(X) = 18 - 0.1 \cdot x\) | \(res.sd = 2.5\), \(res.sd = 2\), \(res.sd = 1.5\) |
| 13 | \(f(X) = 25 + 0.02 \cdot x -0.003 \cdot (x - 45)^2 + 0.00006 \cdot (x - 54)^3\) | \(n = 600\), \(n = 550\), \(n = 500\) |
| 14 | \(f(X) = 7 + 3.5 \cdot \cos{\Big ( {x - 30 \over 12 } \Big )}\) | \(train.percent = 0.9\), \(train.percent = 0.85\), \(train.percent = 0.8\) |
| 15 | \(f(X) = 16 - 0.1 \cdot x\) | \(res.sd = 2\), \(res.sd = 3\), \(res.sd = 4\) |
| 16 | \(f(X) = 25 + 0.02 \cdot x -0.003 \cdot (x - 45)^2 + 0.00006 \cdot (x - 54)^3\) | \(train.percent = 0.3\), \(train.percent = 0.25\), \(train.percent = 0.2\) |
| 17 | \(f(X) = 9 + 3.5 \cdot \cos{\Big ( {x - 30 \over 12 } \Big )}\) | \(n = 450\), \(n = 400\), \(n = 350\) |
| 18 | \(f(X) = 12 + 3.5 \cdot \sin{\Big ( {x - 30 \over 9 } \Big )}\) | \(res.sd = 6\), \(res.sd = 7\), \(res.sd = 8\) |
| 19 | \(f(X) = 15 + 0.02 \cdot x -0.005 \cdot (x - 45)^2 + 0.00006 \cdot (x - 54)^3\) | \(res.sd = 6\), \(res.sd = 7\), \(res.sd = 8\) |
| 20 | \(f(X) = 14 - 0.05 \cdot x\) | \(n = 850\), \(n = 900\), \(n = 950\) |
Источники
- James G., Witten D., Hastie T. and Tibshirani R. An Introduction to Statistical Learning with Applications in R. URL: http://www-bcf.usc.edu/~gareth/ISL/ISLR%20First%20Printing.pdf