Skip to content

Лабораторные работы по предмету "Системы и методы принятия решений"

License

Notifications You must be signed in to change notification settings

IHappyPlant/SMPR

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 
 
 
 
 
 
 
 
 
 
 

Repository files navigation

Теория машинного обучения

Навигация

Метрические алгоритмы классификации

Метрические алгоритмы классификации - алгоритмы, основанные на вычислении оценок сходства между объектами. Одним из вариантов сходства является расстояние между объектами. Тогда вводится функция расстояния .

Сводная таблица по метрическим алгоритмам

Алгоритм Параметры LOO Кол-во ошибок
STOLP-kwNN k = 20, q = 0.1; после отбора эталонов k = 8, q = 0.1 STOLP сам по себе вместо LOO 3
kNN k = 6 0.03 5
kwNN k = 6, q = 1 0.03 5
Потенциальные функции Ядро - квартическое, h(setosa) = 1, h(!setosa) = 0.4, Gamma_i - случайные вместо LOO - подбор Gamma_i 5
Парзеновское окно Ядро - Прямоугольное, h = 0.4 0.04 6
Парзеновское окно Ядро - Треугольное, h = 0.4 0.04 6
Парзеновское окно Ядро - Епанечникова, h = 0.4 0.04 6
Парзеновское окно Ядро - Квартическое, h = 0.4 0.04 6
Парзеновское окно Ядро - Гауссовское, h = 0.1 0.04 6

Метод k ближайших соседей (kNN)

Один из самых простых метрических алгоритмов классификации. Работает следующим образом: пусть дан классифицируемый объект z и обучающая выборка . Требуется определить класс объекта z на основе данных из обучающей выборки. Для этого:

  1. Вся выборка сортируется по возрастанию расстояния от объекта z до каждого объекта выборки.
  2. Проверяются классы k ближайших соседей объекта z. Класс, встречаемый наиболее часто, присваивается объекту z.

Вход алгоритма - классифицируемый объект, обучающая выборка и параметр k - число рассматриваемых ближайших соседей. Выход алгоритма - класс классифицируемого объекта.

При k = 1 алгоритм превращается в 1NN, то есть объекту z присваивается класс его первого ближайшего соседа, все остальные объекты выборки не учитываются.

Оптимальное значение k подбирается по Критерию Скользящего Контроля (LOO). Суть критерия: пусть дана обучающая выборка xl. Требуется определить оптимальное значение параметра k для данной выборки. Для этого:

  1. Из выборки удаляется i-й объект .
  2. Выбранный алгоритм классификации запускается для x_i и оставшейся выборки. По окончании работы алгоритма полученный класс объекта сравнивается с его реальным классом. При их несовпадении сумма ошибки увеличивается на 1.
  3. Шаги 1 и 2 повторяются для каждого объекта выборки при фиксированном k. По окончании работы алгоритма полученная сумма ошибки sum делится на размер выборки l: sum=sum/l . Потом значение k меняется, и алгоритм повторяется для нового значения. k с наименьшим значением суммы ошибки будет оптимальным.

Реализация

При реализации алгоритма, в качестве обучающей выборки использовалась выборка ирисов Фишера. В качестве признаков объектов использовались значения длины и ширины лепестка. Значение k подбиралось по LOO.

Алгоритм:

kNN <- function(xl, z, k) names(which.max(table(sort_objects_by_dist(xl, z)[1:k, ncol(xl)])))

где xl - обучающая выборка.

kNN.png

Достоинства алгоритма:

  1. Простота реализации
  2. Хорошее качество, при правильно подобранной метрике и параметре k

Недостатки алгоритма:

  1. Необходимость хранить выборку целиком
  2. Малый набор параметров
  3. Качество классификации сильно зависит от выбранной метрики

Метод k взвешенных ближайших соседей (kwNN)

Отличается от kNN тем, что вес ближайших соседей зависит не от ранга соседа, а от расстояния до объекта z. В kNN каждый из k ближайших соседей имеет вес равный единице, а все остальные объекты выборки имеют вес, равный нулю. Поэтому можно было говорить о частоте появления класса среди k ближайших соседей.
В методе kwNN для задания весов k ближайшим соседям должна использоваться невозрастающая функция.
Мы будем использовать функцию весов w(i)=q^i , i - ранг соседа, q - подбираемый параметр. Значение q подбирается для каждого k по LOO.

Реализация

При реализации использовалась та же выборка ирисов Фишера. Значения k и q подбирались по LOO
Алгоритм:

w.kwnn <- function(i, k, q) (i <= k) * q^i # Весовая функция
kwNN <- function(xl, z, k, q) {
  ordered_xl <- sort_objects_by_dist(xl, z)
  weights <- w.kwnn(1:nrow(ordered_xl), k, q)
  names(weights) <- ordered_xl[, ncol(ordered_xl)]
  sum_by_class <- sapply(unique(sort(names(weights))), function(class, weights) sum(weights[names(weights) == class]), weights)
  names(which.max(sum_by_class))
}

где xl - обучающая выборка.

kwNN.png

Достоинства алгоритма

  1. Простота реализации
  2. Хорошая эффективность при правильно подобранной метрике, k и q. В общем случае, более высокая эффективность, чем у kNN.

Недостатки. Вообще говоря, те же, что и у kNN:

  1. Необходимость хранить выборку целиком
  2. Малый набор параметров
  3. Сильная зависимость качества классификации от выбранной метрики (особенность всех метрических алгоритмов)

На выборке ирисов Фишера kwNN показывает такую же эффективность, как и kNN. Поэтому, чтобы показать более высокую эффективность kwNN, подберём выборку вруную так, чтобы kNN ошибся при классификации, а kwNN - нет. Приведём пример:

kNN_vs_kwNN.png

При k = 5 точку с координатами (5; 5) kNN отнёс с синему классу, так как синих объектов среди k ближайших соседей больше, а kwNN при k = 5, q = 0.1, отнёс её к классу красных объектов, так как они расположены ближе к точке, и влияют на классификацию сильнее.

Метод парзеновского окна

Идея этого алгоритма в том, чтобы присваивать вес каждому объекту выборки не на основе ранга близости i-го объекта к классифицируемому, а на основе расстояния от классифицируемого объекта до данного объекта выборки. В отличии от метода ближайших соседей, сортировка объектов выборки по расстоянию здесь не нужна.

Функция веса выглядит следующим образом:
Где i - номер объекта выборки, z - классифицируемый объект, xi - i-й объект выборки, h - ширина окна, K - функция ядра.

Функция ядра - произвольная чётная функция, невозрастающая на [0, +inf). На практике применяются следующие функции ядра:

  1. - Прямоугольное ядро
  2. - Треугольное ядро
  3. - Ядро Епанечникова
  4. - Квартическое ядро
  5. - Гауссовское ядро

Где

Программная реализация алгоритма:

parzen <- function(xl, h, distances, kernelFunction = kernel.G) {
  # xl - выборка
  # h - ширина окна
  # distances - расстояния от объекта z до каждого объекта из xl 
  # Расстояния считаются заранее, поэтому сам объект z здесь не нужен
  # kernelFunction - используемая функция ядра. По умолчанию Гауссовское
  l <- nrow(xl)
  n <- ncol(xl)
  classes <- xl[1:l, n] # Классы объектов выборки
  weights <- table(classes) # Таблица для весов классов
  weights[1:length(weights)] <- 0
  for (i in 1:l) { # Для каждого объекта выборки
    class <- xl[i, n] # Берём его класс
    r <- distances[i] / h
    weights[class] <- weights[class] + kernelFunction(r) # И прибавляем его вес к общему весу его класса
  }
  if (max(weights) != 0) # Если точке присвоились какие-нибудь веса классов (точка попала в окно)
    return (names(which.max(weights))) # Вернуть класс с максимальным весом
  return (0) # Иначе вернуть 0
}

Значение h подбиралось в пределах от 0.1 до 2 по LOO. Алгоритм тестировался на выборке ирисов Фишера для разных функций ядра:

Прямоугольное ядро: Треугольное ядро: Ядро Епанечникова: Квартическое ядро: Гауссовское ядро:

У всех ядер, кроме Гауссовского один и тот же недостаток: они не способны проклассифицировать точки, не попавшие в окно (те, для которых |r| > 1).
Значение h имеет смысл подбирать в пределах небольших значений, т. к. при большом значении h слишком много точек попадает в окно, и качество классификации падает крайне значительно.
Для примера приведём карту классификации для Гауссовского ядра при попытке выбора h по LOO в пределах от 1 до l (длина выборки):

Достоинства алгоритма:

  1. Простота реализации
  2. Высокое качество классификации при правильно подобранном h
  3. Для классификации объекта не требуется сортировка выборки, что сокращает время классификации

Недостатки алгоритма:

  1. Необходимость хранить выборку целиком
  2. Малый набор параметров
  3. Если ни один объект выборки не попал в радиус окна h вокруг объекта z, то алгоритм не способен его проклассифицировать. Исправляется использованием Гауссовского ядра и/или переменной ширины окна.

Метод потенциальных функций

Метод получается, если поместить центр парзеновского окна не в классифицируемый объект, а в каждый объект обучающей выборки.
Пусть дана обучающая выборка xl и объект z, который требуется классифицировать. Тогда:

  1. Для каждого объекта выборки задаётся ширина окна h_i (выбирается из собственных соображений).
  2. Для каждого объекта выборки задаётся сила потенциала gamma_i. Алгоритм подбора gamma_i указан далее.
  3. Каждому объекту выборки присваивается вес по формуле , - функция ядра.
  4. Суммируются веса объектов одинаковых классов. Класс с самым большим весом присваивается объекту z.

Алгоритм подбора gamma_i:

  1. По умолчанию все gamma_i задаются равными нулю. Задается максимально допустимое число ошибок.
  2. Из обучающей выборки случайным образом, или последовательно выбирается объект x_i.
  3. Для x_i запускается алгоритм классификации.
  4. Если полученный класс не совпал с реальным, то сила потенциала для выбранного объекта увеличивается на 1. Иначе шаги 2-4 повторяются.
  5. Алгоритм классификации с полученными значениями потенциалов запускается для каждого объекта выборки. Подсчитывается число ошибок.
  6. Если число ошибок меньше заданного на шаге 1, то алгоритм завершает работу. Иначе повторяются шаги 2-6.

Программная реализация метода потенциальных функций:

pF <- function(distances, potentials, h, xl, kernelFunction = kernel.G) {
  l <- nrow(xl)
  n <- ncol(xl)
  classes <- xl[, n]
  weights <- table(classes) # Таблица для весов классов
  weights[1:length(weights)] <- 0 # По умолчанию все веса равны нулю
  for (i in 1:l) { # Для каждого объекта выборки
    class <- xl[i, n] # Берется его класс
    r <- distances[i] / h[i]
    weights[class] <- weights[class] + potentials[i] * kernelFunction(r) # Считается его вес, и прибавляется к общему ввесу его класса
  }
  if (max(weights) != 0) return (names(which.max(weights))) # Если есть веса больше нуля, то вернуть класс с наибольшим весом
  return ("") # Если точка не проклассифицировалась, то вернуть пустую строку
}

Программная реализация алгоритма подбора gamma_i:

getPotentials <- function(xl, h, eps, kernelFunction = kernel.G) {
  # Получить потенциалы всех объектов выборки
  l <- nrow(xl)
  n <- ncol(xl)
  potentials <- rep(0, l)
  err <- eps + 1
  # Посчитаем расстояния от каждого объекта выборки до остальных
  distances <- matrix(0, l, l)
  for (i in 1:l)
    distances[i,] <- getDistances(xl, c(xl[i, 1], xl[i, 2]))
  # Пока число ошибок больше заданного
  while (err > eps) {
    while (TRUE) {
      # Пока не получим несоответствие классов, чтобы обновить потенциалы
      rand <- sample(1:l, 1)
      class <- pF(distances[rand,], potentials, h, xl)
      if (class != xl[rand, n]) {
        potentials[rand] = potentials[rand] + 1
        break
      }
    }
    # Подсчет числа ошибок
    err <- 0
    for (i in 1:l) {
      class <- pF(distances[i,], potentials, h, xl)
      err <- err + (class != xl[i, n])
    }
  }
  return (potentials)
}

Алгоритм тестировался на выборке ирисов Фишера. Использовалось квартическое ядро. Радиусы потенциалов h для класса setosa были положены равными 1, так как класс далеко от остальных, а радиусы для versicolor и virginica 0.4, так как это число было оптимальным значением h для метода парзеновского окна.
Результат работы:

Достоинства алгоритма:

  1. Большое количество параметров для подбора.
  2. После настройки силы потенциалов, объекты выборки с нулевыми потенциалами можно не использовать при классификации. Благодаря этому - высокая скорость работы.

Недостатки алгоритма:

  1. Слишком грубая настройка параметров h_i и gamma_i.
  2. Неопределённое время работы алгоритма подбора gamma_i.
  3. При случайном выборе объектов из выборки при подборе gamma_i, результат работы на одной и той же выборке будет разным.
  4. Исходя из предыдущих пунктов - неопределённое качество классификации.

Алгоритм STOLP

Введём понятие отступа. Отступ - величина, показывающая, насколько объект является типичным представителем своего класса. Отступ равен разности между степенью близости объекта к своему классу и суммой близостей объекта ко всем остальным классам.

Все объекты обучающей выборки можно разделить на несколько типов:

  1. Эталонные объекты - наиболее типичные представители своего класса
  2. Неинформативные - объекты, не влияющие значительным образом на качество классификации
  3. Пограничные - объекты, имеющие отступ, близкий к нулю. Незначительное изменение в выборке, алгоритме, метрике и т.п. может повлиять на их классификацию.
  4. Ошибочные - объекты с отрицательными отступами, классифицируемые неверно.
  5. Шумовые объекты (выбросы) - малая группа объектов с большими отрицательными отступами. Их удаление улучшает качество классификации.

Идея алгоритма STOLP состоит в том, чтобы уменьшить исходную выборку, выбрав из неё эталонные объекты. Такой подход уменьшит размер выборки, и может улучшить качество классификации. На вход подаётся выборка, допустимый порог ошибок и порог фильтрации выбросов. Алгоритм:

  1. Удалить из выборки все выбросы (объекты, отступ которых меньше порога фильтрации выбросов).
  2. Пересчитать все отступы заново, взять по одному эталону из каждого класса (объекты с наибольшим положительным отступом), и добавить их в множество эталонов.
  3. Проклассифицировать объекты обучающей выборки, взяв в качестве обучающей выборки для этого множество эталонов. Посчитать число ошибок.
  4. Если число ошибок меньше заданного порога, то алгоритм завершается.
  5. Иначе присоединить ко множеству эталонов объекты с наименьшим отступом из каждого класса из числа классифицированных неправильно.
  6. Повторять шаги 3-5 до тех пор, пока множество эталонов и обучающая выборка не совпадут, или не сработает проверка в пункте 4.

Протестируем алгоритм на выборке ирисов Фишера.
В качестве алгоритма классификации будем использовать kwNN(k = 20, q = 0.1)
Реализация функции отступа на основе kwNN:

margin <- function(xl, z, k, q, target_class) {
  ordered_xl <- sort_objects_by_dist(xl, z) # Сортировка выборки
  weights <- w.kwnn(1:nrow(ordered_xl), k, q) # Веса объектов выборки
  names(weights) <- ordered_xl[, ncol(ordered_xl)]
  sum(weights[names(weights) == target_class]) - sum(weights[names(weights) != target_class]) # Отступ
}

Результат работы STOLP:
На оригинальной выборке kwNN с такими параметрами отработал с 7 ошибками, а на выборке из эталонов - всего 3. Выборка эталонов составила 8 объектов, значит, выборка уменьшилась 18.75 раз.

Достоинства алгоритма:

  1. Значительно сокращает размер выборки
  2. Улучшает качество классификации
  3. Ускоряет классификацию

Недостатки алгоритма:

  1. Сложность реализации

Байесовские алгоритмы классификации

Байесовские алгоритмы классификации основаны на идее о том, что обучающая выборка формируется из неизвестного вероятностного распределения, имеющего некоторую плотность.
В основе байесовского подхода лежит теорема о том, что если для классов известны вероятности P и плотности распределения p, то можно записать решающее правило в виде
, где - важность класса y,
и оно будет оптимальным, т. е. обладать наименьшей вероятностью ошибки.
На практике вероятности и плотности распределения для классов обычно неизвестны, и их приходится восстанавливать по выборке. Из за этого решающее правило перестаёт быть оптимальным, так как восстановить эти значения по выборке можно только с некоторой погрешностью.
Существуют параметрический и не параметрический подходы к задаче восстановления плотности.

При параметрическом подходе плотность распределения выборки полагается известной с точностью до параметра .

Нормальный дискриминантный анализ

Это метод классификации, основанный на параметрическом подходе, специальный случай, при котором плотности распределения всех классов полагаются многомерными нормальными.
Функция плотности многомерного нормального распределения выглядит следующим образом:
, где
- объект, состоящий из n признаков
- мат. ожидание каждого признака
- матрица ковариации признаков. Симметричная, невырожденная, положительно определённая.

Линии уровня нормального распределения

  1. Если признаки некоррелированы, т. е. матрица ковариации диагональна, то линии уровня имеют форму эллипсоидов, параллельных осям координат, вытянутых относительно признака, значение для которого в матрице выше.

  2. Если признаки коррелированы, то есть матрица ковариации не диагональна, то линии уровня имеют форму эллипсоидов, наклонённых относительно осей координат.

    Сама программа, рисующая линии уровня представлена здесь

Наивный байесовский классификатор

Если предположить, что все признаки объектов выборки сформированы независимо, то значение функции плотности для класса y для объекта x можно представить в виде , где - i-й признак объекта x. Это упрощает задачу, так как оценивать несколько одномерных плотностей легче, чем одну многомерную.
Однако, на практике такая ситуация встречается редко, поэтому алгоритм получил название "Наивный байесовский классификатор". Обычно он используется, как эталон при сравнении различных алгоритмов классификации.
Решающее правило принимает вид:

Примеры работы алгоритма:

В приведённых примерах все признаки действительно взяты независимо, поэтому алгоритм показывает хорошее качество классификации. В практических задачах результат был бы намного хуже.

Достоинства алгоритма:

  1. Простота реализации.
  2. Низкие вычислительные затраты при обучении и классификации.
  3. Если признаки (почти) независимы, то алгоритм (почти) оптимален.

Недостатки алгоритма:

  1. В общем случае - низкое качество классификации.

Приложение, реализующее алгоритм, представлено здесь.

Подстановочный (Plug-in) алгоритм

Алгоритм относится к нормальному дискриминантному анализу и применяется для многомерных (в данном случае - для многомерного нормального) распределений.
Решающее правило имеет вид: , где - восстановленные вероятность и плотность распределения класса соответственно.



где - восстановленная матрица ковариации. Матрица ковариации вычисляется по формуле:

Восстановление мат. ожидания:

get_mu <- function(xm) colMeans(xm)

Восстановление матрицы ковариации:

get_sigma <- function(xm, mu) {
  sum <- 0
  for (i in 1:nrow(xm)) {
    xi <- matrix(c(xm[i,1], xm[i,2]), 1, 2)
    sum <- sum + t(xi - mu) %*% (xi - mu)
  }
  sum / (nrow(xm)-1)
}

Полученные матрицу ковариации и вектор мат. ожиданий подставляют в формулу плотности, которую, в свою очередь, подставляют в решающее правило. Поэтому алгоритм получил название "Подстановочный".

Алгоритм позволяет построить разделяющую кривую между классами, решив уравнение:

Это уравнение, при раскрытии переменных, превратится в уравнение кривой второго порядка, при подстановке точек, в которое, получим разделяющую кривую между классами.

Примеры работы алгоритма:

Приложение, реализующее алгоритм, представлено здесь.

Линейный дискриминант Фишера

ЛДФ является частным случаем подстановочного алгоритма, для случая, когда матрицы ковариации всех классов равны.
В этом случае матрицу ковариации можно восстановить для всех классов сразу, по формуле:
где |Y| - мощность множества классов, x - объекты выборки. Программная реализация:

get_sigma <- function(xm, mu) {
  sum <- 0
  m <- nrow(xm)
  for (i in 1:(m/2)) {
    xi <- matrix(c(xm[i,1], xm[i,2]), 1, 2)
    sum <- sum + t(xi - mu[1,]) %*% (xi - mu[1,])
  }
  for (i in (m/2+1):m) {
    xi <- matrix(c(xm[i,1], xm[i,2]), 1, 2)
    sum <- sum + t(xi - mu[2,]) %*% (xi - mu[2,])
  }
  sum / (m-2)
}

Уравнение разделяющей поверхности упрощается, и превращается из уравнения кривой второго порядка в уравнение прямой.

Примеры работы алгоритма:

Приложение, реализующее алгоритм, представлено здесь.

Сравнение Plug-in и ЛДФ

Очевидно, что подстановочный алгоритм покрывает больше случаев, и является намного более универсальным, чем ЛДФ.
Но ЛДФ создавался специально для работы с особыми случаями, а именно - когда матрицы ковариации классов равны.

Сравним работу Plug-in и ЛДФ на таких данных:
Plug-in:
ЛДФ:
Plug-in:
ЛДФ:
Plug-in:
ЛДФ:

Вывод: для случая, когда матрицы ковариации классов равны, и Plug-in и ЛДФ показывают схожие результаты, что неудивительно, т. к. ЛДФ - частный случай Plug-in. Тем не менее, ЛДФ значительно проще в реализации, и требует меньше вычислительной мощности, благодаря более простой формуле. Также, в случае, когда классы накладываются друг на друга, Plug-in "ломается", а ЛДФ - нет.

EM-алгоритм

Класс не во всех случаях можно описать одним распределением. Иногда класс представляет собой смесь распределений. В этом случае функция правдоподобия p(x) вычисляется по формуле: , где j - номер компоненты смеси, - её апостериорная вероятность. Задача разделения смеси состоит в том, чтобы имея обучающую выборку, выделить из неё компоненты смеси, и оценить вектор параметров
Для разделения смеси применяется EM-алгоритм. Он состоит из двух больших шагов: E (expectation) - шага и M (maximization) - шага. В начале алгоритма задаётся начальное приближение , затем на E - шаге по вычисляется вектор скрытых переменных G. На M - шаге по текущим значениям G и вычисляется новое значение . Процесс повторяется, пока G и не стабилизируются.

E - шаг

На E - шаге вычисляются - апостериорные вероятности того, что объект принадлежит j-й компоненте смеси.

M - шаг

На M - шаге решается задача максимизации правдоподобия:
Путём некоторых сложных преобразований она разбивается на k подзадач (по количеству компонент):

Рассмотрим смесь многомерных гауссовских распределений. Для простоты предположим, что все компоненты имеют диагональные матрицы ковариации. Тогда функция правдоподобия объекта x для класса y вычисляется по формуле:


И решение задачи максимизации правдоподобия на M - шаге имеет следующий вид:

EM - алгоритм с фиксированным числом компонент:

Вход: - обучающая выборка, число компонент смеси, начальное приближение вектора параметров, параметр критерия останова
Выход: - оптимизированный вектор параметров

  1. E - шаг: для всех :
    ;
  2. M - шаг: для всех j = 1,...,k вычислить оптимальные мат. ожидание и дисперсию, пересчитать (формулы выше).
  3. Если (), то перейти на п. 1
  4. Вернуть

Если количество компонент неизвестно, то применяется модификация: EM - алгоритм с последовательным добавлением компонент.

EM - алгоритм с последовательным добавлением компонент

Вход: - обучающая выборка, минимальное число объектов для восстановления компоненты, параметр разброса правдоподобия объектов, параметр критерия останова
Выход: - оптимизированный вектор параметров

  1. Старт: восстановить мат ожидание и дисперсию по всей выборке,
  2. Выделить объекты с низким правдоподобием:
  3. Если , то переход к п. 6
  4. k := k + 1. Добавим новую компоненту, и вычислим начальное приближение для неё: восстановим мат. ожидание и матрицу ковариации для U и пересчитаем
  5. EM
  6. Вернуть

Программная реализация вышеприведённых алгоритмов:

EM_seq <- function(xl, m0, r, delta) {
  # EM - алгоритм с последовательным добавлением компонент
  l <- nrow(xl)
  objects <- xl[,-ncol(xl)]
  mu <- matrix(get_mu(objects), 1, 2)
  sigma <- get_sigma(objects, mu)
  theta <- cbind(sigma, t(mu))
  k <- 1
  w <- 1
  max_iter <- 7 # максимум компонент, чтобы не зацикливался
  for (iter in 1:max_iter) {
    phi_all <- get_phi_all(objects, theta, w)
    max_r <- max(phi_all) * r
    u <- which(phi_all < max_r)
    u_capacity <- length(u)
    
    if (u_capacity < m0) break
    
    k <- k + 1
    wk <- u_capacity / l
    w <- sapply(1:length(w), function(i) w[i] * (1 - wk))
    w <- c(w, wk)
    mu <- matrix(get_mu(objects[u,]), 1, 2)
    sigma <- get_sigma(objects[u,], mu)
    theta <- rbind(theta, cbind(sigma, t(mu)))
    theta <- EM(xl, k, theta, delta, w)
    w <- c(theta[1:k,4]) # Уберём w из полученного theta
    theta <- theta[,-4]
  }
  return (cbind(theta, as.matrix(c(w, rep(0, k)))))
}

EM <- function(xl, k, theta, delta, w) {
  # EM - алгоритм с фиксированным числом компонент
  l <- nrow(xl)
  objects <- xl[,-ncol(xl)]
  n <- ncol(objects)
  g <- matrix(0, l, k)
  g0 <- matrix(0, l, k)
  max_iter <- 100 # максимум итераций, чтобы не зацикливался
  for (iter in 1:max_iter) {
    # E - шаг
    for (i in 1:l)
      for (j in 1:k) {
        g0[i, j] <- g[i, j]
        tmp <- sum(sapply(1:k, function(s) w[s] * phi(objects[i,], theta[(2*s-1):(2*s), 3], theta[(2*s-1):(2*s), (1:2)])))
        tmp1 <- w[j] * phi(objects[i,], theta[(2*j-1):(2*j), 3], theta[(2*j-1):(2*j), (1:2)])
        g[i, j] <- tmp1 / tmp 
      }
    # Жуткий костыль: если в какой-то ячейке nan (почему то иногда плотность становится нулевой и tmp становится нулём), тогда запишем в неё результат из g0
    g <- check_fix_nan(g0, g)
    # M - шаг
    for (j in 1:k) {
      w[j] <- sum(g[,j]) / l
      mu <- sapply(1:n, function(i) sum(g[,j] * objects[,i]) / (l * w[j])) # Оптимальное мат ожидание
      sigma <- matrix(0, 2, 2)
      sigma[1,1] <- sum(g[,j] * (objects[,1] - mu[1])^2) / (l * w[j]) # Оптимальная матрица ковариации (диагональная)
      sigma[2,2] <- sum(g[,j] * (objects[,2] - mu[2])^2) / (l * w[j])
      theta[(2*j-1),] <- c(sigma[1,], mu[1])
      theta[(2*j),] <- c(sigma[2,], mu[2])
    }
    if (to_stop(g, g0) <= delta) break
  }
  return (cbind(theta, as.matrix(c(w, rep(0, k)))))
}

Для построения разделяющей поверхности между классами будем использовать сеть радиальных базисных функций. Пусть классы представляют собой смесь нормальных распределений с диагональными матрицами ковариации. Тогда функция правдоподобия записывается в виде: , где N - нормировочный множитель,
- взвешенная евклидова метрика
В данном случае плотность распределения означает расстояние от объекта x до центра j-й компоненты (до ). Т. к. функции зависят от расстояния между x и фиксированными точками пространства, их называют радиальными. Алгоритм классификации имеет вид:
Алгоритм состоит из трёх этапов: на первом вычисляются функции правдоподобия для всех компонент всех классов. На втором происходит умножение каждой функции правдоподобия на вероятность своей компоненты и суммирование по классам. На третьем этапе происходит домножение полученных сумм на соответствующие и , сравнение полученных выражений и возврат
Примеры работы алгоритма:




Приложение, реализующее алгоритм, представлено здесь

Линейные алгоритмы классификации

Рассмотрим задачу классификации с множеством объектов и множеством ответов Y = {-1, +1}.
Рассмотрим модель алгоритма , где w - вектор параметров, а f(x, w) - дискриминантная функция.
Если f(x, w) > 0, то алгоритм отнесёт объект к классу +1, иначе - к классу -1.
Уравнение f(x, w) = 0 задаёт разделяющую поверхность.
Если , то алгоритм классификации называется линейным, и решающее правило принимает вид:

Обучение состоит в минимизации аппроксимированного эмпирического риска:
, где - функция потерь, M - отступ
и заключается в подборе компонент вектора w.

Для минимизации Q(w) применим модифицированный метод градиентного спуска.

Метод стохастического градиента

Суть метода: выбирается начальное приближение вектора w, и запускается итерационный процесс, при котором на каждом шаге вектор w изменяется в направлении наиболее быстрого убывания Q. Это направление противоположно направлению вектора градиента : , где - темп обучения. В этом методе прецеденты (пары "объект-ответ") для пересчёта w выбираются случайно, поэтому метод получил название "Метод стохастического градиента".
Алгоритм:

  1. Инициализировать вектор w.
  2. Вычислить начальное значение функционала .
  3. Выбрать из выборки объект .
  4. Вычислить ошибку алгоритма на выбранном объекте: .
  5. Пересчитать значения .
  6. Пересчитать значение .
  7. Повторять п.3 - п.6 до тех пор, пока значение Q и/или вектор w не перестанут изменятся.

Программная реализация:

gradient <- function(xl, eta, lambda, rule, loss_function, eps_q, stop_by_margins) {
  l <- nrow(xl)
  n <- ncol(xl)
  w <- matrix(c(runif(n-1, -1/(2*(n-1)), 1/(2*(n-1)))), 1, 3) # Инициализация w случайными значениями
  objects <- xl[,-n]
  classes <- xl[, n]
  q <- sum(sapply(1:l, function(i) loss_function(margin(w, objects[i,], classes[i])))) # Начальная инициализация q
  q_full <- matrix(q, 1, 1)
  cnt <- 0
  while (T) {
    cnt <- cnt + 1 # Счётчик итераций
    margins <- sapply(1:l, function(i) margin(w[cnt,], objects[i,], classes[i]))
    errors <- which(margins < 0)
    
    if (length(errors) == 0 && stop_by_margins == T) break;
    
    if (length(errors) > 0) rand <- sample(errors, 1)
    else rand <- sample(1:l, 1)
    
    eps <- loss_function(margin(w[cnt,], objects[rand,], classes[rand])) # Ошибка алгоритма на объекте
    
    eta <- 1 / (objects[rand,] %*% objects[rand,])^2 # Пересчёт темпа обучения
    if (length(errors) == 0) eta <- eta / 2
    
    # Обновление весов
    w <- rbind(w, rule(w[cnt,], objects[rand,], classes[rand], eta))
    
    # Пересчёт q
    q_prev <- q
    q <- (1 - lambda) * q + lambda * eps
    q_full <- rbind(q_full, q)
    
    if (abs(q_prev - q) / max(q_prev, q) <= eps_q) { print("exit by q"); break; }
    else if (cnt == 20000) { print("exit by cnt"); break; }
    
  }
  w <- cbind(w, q_full)
  w
}

Преимущества метода:

  1. Простота реализации.
  2. Метод легко обобщается и на нелинейные классификаторы.
  3. Позволяет быстро настраивать веса даже на избыточно больших выборках.

Недостатки:

  1. Неустойчивая сходимость (может сойтись к локальному экстремуму, может сходится очень медленно, а может вообще не сходится).
  2. При большой размерности пространства или малой длине выборки возможно переобучение. Тогда классификация становится неустойчивой.

Метод ADALINE

Возьмём . Вычислив производную по w, получим: .
Получим правило обновления весов на каждой итерации стохастического градиента: .
Это правило называется дельта-правилом и реализуется программно следующим образом:

adaline.get_w <- function(w, object, class, eta)  w - c(eta) * (w %*% object - class) %*% object

Примеры работы классификатора, обученного методом стохастического градиента с применением дельта-правила:


Приложение, реализующее ADALINE, представлено здесь.

Правило Хебба (персептрон Розенблатта)

Этот принцип обучения классификатора основан на принципах нейрофизиологии. Предполагается, что если нейрон угадывает правильный ответ, то его синапсы усиливаются, а если он часто ошибается, или вообще не используется, то синапсы начинают затухать.
В приложении к этого к классификатору роль синапсов играет вектор w, и правило его обновления формализуется следующим образом:
Функция потерь полагается кусочно-линейной:
Примеры работы персептрона:


Приложение, реализующее метод, представлено здесь.

Логистическая регрессия

Является одновременно линейным и байесовским классификатором. Использует логарифмическую (логистическую) функцию потерь: .
Алгоритм назван так потому, что обычно для его обучения применяется не метод стохастического градиента, а метод Ньютона-Рапсона, в котором решается несколько задач восстановления регрессии. Правило обновления весов: , где - сигмоидная функция. С вероятностной точки зрения, вычисляет апостериорную веростность .
Пример классификации методом логистической регрессии:


Достоинства метода:

  1. Как правило, даёт более качественный результат по сравнению с ЛДФ, т. к. основано на менее жестких предположениях, а так же по сравнению с ADALINE и правилом Хебба, т. к. использует более плавную функцию потерь.
  2. Возможность оценивания апостериорной вероятности.

Недостатки метода:

  1. Если не выполняются предположения теоремы, то качество работы может ухудшиться.

Приложение, реализующее метод логистической регрессии, представлено здесь

Сравнение линейных классификаторов

Сравним работу всех трёх вышеприведённых классификаторов. Не будем иллюстрировать процесс схождения, чтобы не перегружать картинку, покажем только конечный результат:


Вывод: ADALINE и Логистическая регрессия выигрывают в точности у правила Хебба, когда выборка линейно разделима. Но когда выборка плохо разделима, или не разделима, ADALINE проигрывает Логистической регрессии и правилу Хебба. Следовательно, алгоритм классификации для конкретной задачи в зависимости от исходных данных.

Support Vector Machine (метод опорных векторов)

SVM - алгоритм классификации, считающийся одним из самых лучших. Он назван так, потому что положение разделяющей классы гиперплоскости зависит лишь от небольшого количества элементов обучающей выборки. Они и называются опорными векторами.
С помощью функции ядра метод обобщается на случай нелинейных разделяющих поверхностей. По сути вид поверхности зависит от используемой функции ядра.
Если выборка линейно разделима, то почти всегда провести разделяющую поверхность можно не единственным образом. В SVM она выбирается таким образом, чтобы она отстояла максимально далеко от обоих классов. Умножим алгоритм a(x) на некоторую константу так, чтобы минимальный отступ в каждом классе был равен единице. Тогда полосу, разделяющую классы описывает множество точек: .
Границами полосы служат две гиперплоскости, на которых лежат объекты с минимальным отступом. Ширина полосы вычисляется по формуле:
Для линейно неразделимой выборки разрешим алгоритму допускать некоторое число ошибок при классификации (т. е. разрешим объектам попадать внутрь разделяющей полосы). Для этого введём дополнительные переменные , означающие величину ошибки на объекте i. Тогда выражение (мы хотим минимизировать норму вектора w) превратится в , где C - управляющий параметр, позволяющий находить компромисс между максимизацией ширины полосы и минимизацией суммарной ошибки.

Для реализации алгоритма SVM будет использовать библиотеку kernlab. Из неё мы будем использовать функцию ksvm. У неё очень много различных аргументов, мы будем передавать туда обучающую выборку, разделённую на объекты и ответы, в качестве функции ядра будем использовать "vanilladot", т. к. это ядро отвечает за линейную разделяющую поверхность, и управляющий параметр C. Приведём примеры работы алгоритма для разных случаев, и покажем, как значение C влияет на классификацию:

  1. Если выборка хорошо разделима, то значение C ни на что не влияет (ошибок всё равно нет).

  2. Если классы находятся близко друг к другу, и значение C выставлено большим, то объекты могут попадать внутрь разделяющей полосы.

  3. Если выборка линейно неразделима, то C влияет на количество объектов внутри полосы. При больших значениях C, внутри полосы могут оказываться объекты, не являющиеся опорными.

Достоинства алгоритма:

  1. Высокая скорость работы.
  2. Хорошее качество классификации.

Недостатки:

  1. Сложность реализации (если самому писать, а не использовать готовый набор).
  2. Не существует однозначных рекомендаций для выбора ядра для конкретной задачи.

Приложение, демонстрирующее работу SVM доступно по ссылке.

В начало