R Notebook

Математическое моделирование

Практика 2

Оценка точности модели с дискретной зависимой переменной (Y)

В практических примерах ниже показано:

• как рассчитать матрицу неточностей
• как считать показатели качества модели по матрице неточностей
• как пользоваться наивным байесовским классификатором
• как пользоваться методом kNN (k ближайших соседей)

Модели: наивный байесовский классификатор, kNN (метод k ближайших соседей).
Данные: сгенерированные.

Нам понадобится несколько пакетов для работы с перечисленными методами классификации.

library('class')          # функция knn()
library('e1071')          # функция naiveBayes()
library('MASS')           # функция mvrnorm()
library('emdbook')        # функция dmvnorm()

# ядро
my.seed <- 12345


# Генерируем данные ------------------------------------------------------------

#  Пример 2 ....................................................................
n <- 100               # наблюдений всего
train.percent <- 0.85  # доля обучающей выборки

# фактические значения объясняющих переменных (нормальный закон)
set.seed(my.seed)
x1 <- rnorm(20, 3.7, n = n)

set.seed(my.seed + 1)
x2 <- rnorm(50, 3.3, n = n)

# истинные дискриминирующие правила
rules <- function(x1, x2){
    ifelse((x1 > 20 & x2 < 50) | (x1 < 18 & x2 > 52), 1, 0)
}
# Конец данных примера 2 .......................................................

Пример из лекции: классы не разделяются прямой границей

Рассмотрим пример 2 из лекции:

• \(n = 100\), доля обучающей выборки: 85%
  
• \(\mathrm{x}_1 \sim N(20, 3.7^2)\)
  
• \(\mathrm{x}_2 \sim N(50, 3.3^2)\)
  
• истинные классифицирующие правила:

\[y_i = \begin{cases} 1, & x_{i1} > 20 \, \mathrm{и} \, x_{i2} < 50, \\ 0, & x_{i1} < 18 \, \mathrm{и} \, x_{i2} > 52 \end{cases}\]

# Отбираем наблюдения в обучающую выборку --------------------------------------
set.seed(my.seed)
inTrain <- sample(seq_along(x1), train.percent*n)
x1.train <- x1[inTrain]
x2.train <- x2[inTrain]
x1.test <- x1[-inTrain]
x2.test <- x2[-inTrain]

# используем истинные правила, чтобы присвоить фактические классы
y.train <- rules(x1.train, x2.train)
y.test <- rules(x1.test, x2.test)

# фрейм с обучающей выборкой
df.train.1 <- data.frame(x1 = x1.train, x2 = x2.train, y = y.train)
# фрейм с тестовой выборкой
df.test.1 <- data.frame(x1 = x1.test, x2 = x2.test)

Нарисуем обучающую выборку на графике. Сеткой точек показаны области классов, соответствующие истинным дискриминирующим правилам.

# Рисуем обучающую выборку графике ---------------------------------------------

# для сетки (истинных областей классов): целочисленные значения x1, x2
x1.grid <- rep(seq(floor(min(x1)), ceiling(max(x1)), by = 1),
               ceiling(max(x2)) - floor(min(x2)) + 1)
x2.grid <- rep(seq(floor(min(x2)), ceiling(max(x2)), by = 1),
               each = ceiling(max(x1)) - floor(min(x1)) + 1)

# классы для наблюдений сетки
y.grid <- rules(x1.grid, x2.grid)

# фрейм для сетки
df.grid.1 <- data.frame(x1 = x1.grid, x2 = x2.grid, y = y.grid)

# цвета для графиков
cls <- c('blue', 'orange')
cls.t <- c(rgb(0, 0, 1, alpha = 0.5), rgb(1, 0.5, 0, alpha = 0.5))

# график истинных классов
plot(df.grid.1$x1, df.grid.1$x2, 
     pch = '·', col = cls[df.grid.1[, 'y'] + 1],
     xlab = 'X1', ylab = 'Y1',
     main = 'Обучающая выборка, факт')
# точки фактических наблюдений
points(df.train.1$x1, df.train.1$x2,
       pch = 21, bg = cls.t[df.train.1[, 'y'] + 1], 
       col = cls.t[df.train.1[, 'y'] + 1])

Обучим модель наивного байесовского классификатора и оценим её точность (верность) на обучающей выборке.

# Байесовский классификатор ----------------------------------------------------
#  наивный байес: непрерывные объясняющие переменные

# строим модель
nb <- naiveBayes(y ~ ., data = df.train.1)
# получаем модельные значения на обучающей выборке как классы
y.nb.train <- ifelse(predict(nb, df.train.1[, -3], 
                             type = "raw")[, 2] > 0.5, 1, 0)

# график истинных классов
plot(df.grid.1$x1, df.grid.1$x2, 
     pch = '·',  col = cls[df.grid.1[, 'y'] + 1], 
     xlab = 'X1', ylab = 'Y1',
     main = 'Обучающая выборка, модель naiveBayes')
# точки наблюдений, предсказанных по модели
points(df.train.1$x1, df.train.1$x2, 
       pch = 21, bg = cls.t[y.nb.train + 1], 
       col = cls.t[y.nb.train + 1])

# матрица неточностей на обучающей выборке
tbl <- table(y.train, y.nb.train)
tbl

##        y.nb.train
## y.train  0  1
##       0 51  2
##       1 15 17

# точность, или верность (Accuracy)
Acc <- sum(diag(tbl)) / sum(tbl)
Acc

## [1] 0.8

Как можно видеть на графике, ту часть жёлтого класса, которая расположена в левой верхней области пространства координат, модель классифицирует неверно. Таким образом, байесовская решающая граница не моделирует разрыв жёлтого класса синим. Это происходит потому, что в непрерывном случае наивный байесовский метод исходит из допущения о линейной разделимости двух классов и нормальности распределения объясняющих переменных в них. Однако в этом примере это допущение не выполняется.
Сделаем прогноз классов Y на тестовую выборку и оценим точность модели. Как можно убедиться, точность на тестовой оказывается ниже, чем на обучающей выборке. Учитывая, как ведёт себя классификатор на обучающей выборке, такой модели доверять не стоит.

# прогноз на тестовую выборку
y.nb.test <- ifelse(predict(nb, df.test.1, type = "raw")[, 2] > 0.5, 1, 0)

# матрица неточностей на тестовой выборке
tbl <- table(y.test, y.nb.test)
tbl

##       y.nb.test
## y.test 0 1
##      0 9 0
##      1 3 3

# точность, или верность (Accuracy)
Acc <- sum(diag(tbl)) / sum(tbl)
Acc

## [1] 0.8

Построим модель kNN. Это “ленивый” классификатор, ему не требуется предварительное обучение. А ещё это непараметрический метод, и чем меньше количество ближайших соседей \(k\), тем гибче ведёт себя разделяющая граница. Метод хорошо работает с линейно неразделимыми классами. Зададим \(k = 3\).

# Метод kNN --------------------------------------------------------------------
#  k = 3

# строим модель и делаем прогноз
y.knn.train <- knn(train = scale(df.train.1[, -3]), 
                   test = scale(df.train.1[, -3]),
                   cl = df.train.1$y, k = 3)

# график истинных классов
plot(df.grid.1$x1, df.grid.1$x2, 
     pch = '·', col = cls[df.grid.1[, 'y'] + 1],
     xlab = 'X1', ylab = 'Y1',
     main = 'Обучающая выборка, модель kNN')
# точки наблюдений, предсказанных по модели
points(df.train.1$x1, df.train.1$x2, 
       pch = 21, bg = cls.t[as.numeric(y.knn.train)], 
       col = cls.t[as.numeric(y.knn.train)])

# матрица неточностей на обучающей выборке
tbl <- table(y.train, y.knn.train)
tbl

##        y.knn.train
## y.train  0  1
##       0 52  1
##       1  0 32

# точность (Accuracy)
Acc <- sum(diag(tbl)) / sum(tbl)
Acc

## [1] 0.9882353

Можно видеть, что классификация обучающей выборки методом kNN не отличается от фактических классов наблюдений.
Оценим также точность модели на тестовой выборке.

# прогноз на тестовую выборку
y.knn.test <- knn(train = scale(df.train.1[, -3]), 
                 test = scale(df.test.1[, -3]),
                 cl = df.train.1$y, k = 3)

# матрица неточностей на тестовой выборке
tbl <- table(y.test, y.knn.test)
tbl

##       y.knn.test
## y.test 0 1
##      0 9 0
##      1 0 6

# точность (Accuracy)
Acc <- sum(diag(tbl)) / sum(tbl)
Acc

## [1] 1

Модель kNN оказалась весьма точной на этих данных, чего не скажешь о байесовском классификаторе.

Пример из лекции: классы нормальны и линейно разделимы

Пример 3 из лекции.
- \(n = 100\), доля обучающей выборки: 85%
- класс \(Y=0\): \(X \sim N((23, 49), \begin{pmatrix} 3.5^2 & 0 \\ 0 & 3.4^2 \end{pmatrix})\)
- класс \(Y=1\): \(X \sim N((15, 51), \begin{pmatrix} 2^2 & 0 \\ 0 & 2.5^2 \end{pmatrix})\)

# Генерируем данные ------------------------------------------------------------

# Данные примера 3 .............................................................
n <- 100               # наблюдений всего
train.percent <- 0.85  # доля обучающей выборки

# x-ы -- двумерные нормальные случайные величины
set.seed(my.seed)
class.0 <- mvrnorm(45, mu = c(23, 49), 
                   Sigma = matrix(c(3.5^2, 0, 0, 3.4^2), 2, 2, 
                                  byrow = T))

set.seed(my.seed + 1)
class.1 <- mvrnorm(55, mu = c(15, 51), 
                   Sigma = matrix(c(2^2, 0, 0, 2.5^2), 2, 2, 
                                  byrow = T))

# записываем x-ы в единые векторы (объединяем классы 0 и 1)
x1 <- c(class.0[, 1], class.1[, 1])
x2 <- c(class.0[, 2], class.1[, 2])

# фактические классы Y
y <- c(rep(0, nrow(class.0)), rep(1, nrow(class.1)))

# классы для наблюдений сетки
rules.mv <- function(v.x, v.mean.y0, v.mean.y1, m.sigma.y0, m.sigma.y1){
    ifelse(dmvnorm(v.x, v.mean.y0, m.sigma.y0) > 
               dmvnorm(v.x, v.mean.y1, m.sigma.y1), 0, 1)
}
# Конец данных примера 3 .......................................................


# Отбираем наблюдения в обучающую выборку --------------------------------------
set.seed(my.seed)
inTrain <- sample(seq_along(x1), train.percent*n)
x1.train <- x1[inTrain]
x2.train <- x2[inTrain]
x1.test <- x1[-inTrain]
x2.test <- x2[-inTrain]

# используем истинные правила, чтобы присвоить фактические классы
y.train <- y[inTrain]
y.test <- y[-inTrain]

# фрейм с обучающей выборкой
df.train.1 <- data.frame(x1 = x1.train, x2 = x2.train, y = y.train)
# фрейм с тестовой выборкой
df.test.1 <- data.frame(x1 = x1.test, x2 = x2.test)

Нарисуем обучающую выборку на графике. Сеткой точек показаны области классов, соответствующие истинным дискриминирующим правилам. Это правило создаём, зная истинные законы распределения классов, как максимум из двух плотностей распределения (плотность многомерного закона считаем функцией dmvnorm(), классы точкам сетки присваиваем пользовательской функцией rules.mv()).

# Рисуем обучающую выборку графике ---------------------------------------------

# для сетки (истинных областей классов): целочисленные значения x1, x2
x1.grid <- rep(seq(floor(min(x1)), ceiling(max(x1)), by = 1),
               ceiling(max(x2)) - floor(min(x2)) + 1)
x2.grid <- rep(seq(floor(min(x2)), ceiling(max(x2)), by = 1),
               each = ceiling(max(x1)) - floor(min(x1)) + 1)

# классы для наблюдений сетки
y.grid <- rules.mv(as.matrix(cbind(x1.grid, x2.grid)),
                   c(23, 49), c(15, 51), 
                   matrix(c(3.5^2, 0, 0, 3.4^2), 2, 2, byrow = T),
                   matrix(c(2^2, 0, 0, 2.5^2), 2, 2, byrow = T))

# фрейм для сетки
df.grid.1 <- data.frame(x1 = x1.grid, x2 = x2.grid, y = y.grid)

# цвета для графиков
cls <- c('blue', 'orange')
cls.t <- c(rgb(0, 0, 1, alpha = 0.5), rgb(1,0.5,0, alpha = 0.5))

# график истинных классов
plot(df.grid.1$x1, df.grid.1$x2, 
     pch = '·', col = cls[df.grid.1[, 'y'] + 1],
     xlab = 'X1', ylab = 'Y1',
     main = 'Обучающая выборка, факт')
# точки фактических наблюдений
points(df.train.1$x1, df.train.1$x2,
       pch = 21, bg = cls.t[df.train.1[, 'y'] + 1], 
       col = cls.t[df.train.1[, 'y'] + 1])

Обучим модель наивного байесовского классификатора и оценим её точность (верность) на обучающей выборке. Поскольку объясняющие переменные для классов сгенерированы как двумерные нормальные распределения и сами классы не перекрываются, следует ожидать, что эта модель окажется точной.

# Байесовский классификатор ----------------------------------------------------
#  наивный байес: непрерывные объясняющие переменные

# строим модель
nb <- naiveBayes(y ~ ., data = df.train.1)
# получаем модельные значения на обучающей выборке как классы
y.nb.train <- ifelse(predict(nb, df.train.1[, -3], 
                             type = "raw")[, 2] > 0.5, 1, 0)

# график истинных классов
plot(df.grid.1$x1, df.grid.1$x2, 
     pch = '·',  col = cls[df.grid.1[, 'y'] + 1], 
     xlab = 'X1', ylab = 'Y1',
     main = 'Обучающая выборка, модель naiveBayes')
# точки наблюдений, предсказанных по модели
points(df.train.1$x1, df.train.1$x2, 
       pch = 21, bg = cls.t[y.nb.train + 1], 
       col = cls.t[y.nb.train + 1])

# матрица неточностей на обучающей выборке
tbl <- table(y.train, y.nb.train)
tbl

##        y.nb.train
## y.train  0  1
##       0 34  5
##       1  2 44

# точность, или верность (Accuracy)
Acc <- sum(diag(tbl)) / sum(tbl)
Acc

## [1] 0.9176471

Точность на обучающей выборке высокая. Сделаем прогноз классов Y на тестовую выборку и оценим точность модели.

# прогноз на тестовую выборку
y.nb.test <- ifelse(predict(nb, df.test.1, type = "raw")[, 2] > 0.5, 1, 0)

# матрица неточностей на тестовой выборке
tbl <- table(y.test, y.nb.test)
tbl

##       y.nb.test
## y.test 0 1
##      0 6 0
##      1 1 8

# точность, или верность (Accuracy)
Acc <- sum(diag(tbl)) / sum(tbl)
Acc

## [1] 0.9333333

Наивный байесовский метод разделяет классы на обучающей выборке, ошибаясь лишь в одном наблюдении.

Построим модель kNN. С этими данными у метода не должно возникнуть проблем, так как он не проводит чёткой границы между классами, а в каждом случае ориентируется на соседние наблюдения.

# Метод kNN --------------------------------------------------------------------
#  k = 3

# строим модель и делаем прогноз
y.knn.train <- knn(train = scale(df.train.1[, -3]), 
                   test = scale(df.train.1[, -3]),
                   cl = df.train.1$y, k = 3)

# график истинных классов
plot(df.grid.1$x1, df.grid.1$x2, 
     pch = '·', col = cls[df.grid.1[, 'y'] + 1],
     xlab = 'X1', ylab = 'Y1',
     main = 'Обучающая выборка, модель kNN')
# точки наблюдений, предсказанных по модели
points(df.train.1$x1, df.train.1$x2, 
       pch = 21, bg = cls.t[as.numeric(y.knn.train)], 
       col = cls.t[as.numeric(y.knn.train)])

# матрица неточностей на обучающей выборке
tbl <- table(y.train, y.knn.train)
tbl

##        y.knn.train
## y.train  0  1
##       0 34  5
##       1  2 44

# точность (Accuracy)
Acc <- sum(diag(tbl)) / sum(tbl)
Acc

## [1] 0.9176471

Так и есть. Точность на обучающей выборке близка к идеальной.
Оценка точности на тестовой выборке также показывает, что модель классифицирует верно все наблюдения, кроме двух, близких к границе разделения классов.

# прогноз на тестовую выборку
y.knn.test <- knn(train = scale(df.train.1[, -3]), 
                 test = scale(df.test.1[, -3]),
                 cl = df.train.1$y, k = 3)

# матрица неточностей на тестовой выборке
tbl <- table(y.test, y.knn.test)
tbl

##       y.knn.test
## y.test 0 1
##      0 6 0
##      1 0 9

# точность (Accuracy)
Acc <- sum(diag(tbl)) / sum(tbl)
Acc

## [1] 1

Упражнение 2

1. Построить модели на данных примера 3 с параметрами распределений, соответствующими своему варианту. На графики нанести сетку истинных классов. Определить, какой из методов срабатывает на этих данных лучше, и почему.

2. По матрице неточностей той модели, которая оказалась лучше по \(Acc\), рассчитать характеристики качества и ошибки из лекции: \(TPR\), \(SPC\), \(PPV\), \(NPV\), \(FNR\), \(FPR\), \(FDR\), \(MCC\).

3. Выполненные задачи 1-2 разместить в одном отчёте в репозитории на github.com, выслать ссылку на него на почту преподавателя. В репозитории должны лежать:

• файл README.md с кратким описанием содержимого репозитория;
  
• скрипт генерации отчёта: файл .Rmd в кодировке UTF-8;
  
• графики, сгенерированные в задании.

В отчёте с решением должны присутствовать, кроме блоков кода R, вводный текст с постановкой задачи и выводы (ответ на поставленный вопрос).

Ядро генератора случайных чисел – номер варианта (номер в списке группы). Функцию set.seed() вызывать перед каждой генерацией случайных величин для исходных данных и перед каждой классификации по k средним функцией knn().

Варианты

Все условия, не упомянутые в таблице (величина выборки, закон распределения \(X\) и т.д.) брать из второй практики.

• класс \(Y=0\): \(X \sim N((x1.mean.y0, x2.mean.y0), \begin{pmatrix} x1.sd.y0^2 & 0 \\ 0 & x2.sd.y0^2 \end{pmatrix})\)
  
• класс \(Y=1\): \(X \sim N((x1.mean.y1, x2.mean.y1), \begin{pmatrix} x1.sd.y1^2 & 0 \\ 0 & x2.sd.y1^2 \end{pmatrix})\)

Источники

1. 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