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

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

Как известно, формула Байеса позволяет поменять местами причину и следствие в формулах с условными вероятностями:

\[\mathbb{P}(A \mid B) = \frac{\mathbb{P}(B \mid A)\, \mathbb{P}(A)}{\mathbb{P}(B)}.\]

Пусть в качестве \(A\) у нас будет рассматриваться некоторое множество классов \(C\), которым могут принадлежать наши объекты; а в качестве \(B\) – набор признаков (вектор предикторов) \(X_1, X_2,\ldots, X_n\), по значениям которым мы будем пытаться отнести любой объект к одному из классов. В этих обозначениях формула Байеса перепишется следующим образом:

\[\mathbb{P}(C \mid X_1, X_2,\ldots, X_n) = \frac{\mathbb{P}(X_1, X_2,\ldots, X_n \mid C)\, \mathbb{P}(C)}{\mathbb{P}(X_1, X_2,\ldots, X_n)}.\]

Если, например, рассматриваются три класса: “флет”, “тренд вверх”, “тренд вниз”, а в качестве предикторов (наблюдаемого набора признаков) выступает день недели и вчерашнее состояние тренда, то формула Байеса примет вид:

\[\mathbb{P}(C \mid X_1, X_2) = \frac{\mathbb{P}(X_1, X_2 \mid C)\, \mathbb{P}(C)}{\mathbb{P}(X_1, X_2)},\]

где случайная величина \(C\) может принимать три значения 0 (флет), 1 (uptrend, тренд вверх) и -1 (downtrend, тренд вниз); предиктор \(X_1\) принимает значения Mn, Tu, Wd, Th и Fr (рабочие дни недели с понедельника до пятницы); предиктор \(X_2\) – те же значения, что и \(C\) (флет, тренд вверх или тренд вниз, но имеется в виду направление рынка в течение предыдущего дня).

В числителе \(\mathbb{P}(C)\) – это априорные вероятности того, что любой из рассматриваемых объектов (с пока неизвестными признаками) окажется принадлежащим данному классу. В нашем примере три таких вероятности, для флета и для трендов обоих направлений: \(P_0\), \(P_{1}\), и \(P_{-1}\).

Знаменатель в формуле Байеса \(\mathbb{P}(X_1, X_2)\) – вероятности того, что наугад взятый объект будет иметь заданный набор значений признаков. В нашем примере 15 таких вероятностей:

\[P(\text{Mn}, \text{-1}),\;\; P(\text{Tu}, \text{-1}), \ldots, P(\text{Mn}, \text{0}), \ldots, P(\text{Fr}, \text{1}).\]

\(\mathbb{P}(X_1, X_2 \mid C)\) – это условные вероятности того, что объект, принадлежащий данному классу, имеет заданный набор признаков. В нашем примере \(5\cdot 3\cdot 3 = 45\) таких вероятностей:

\[P(\text{Mn}, -1 \mid 1),\; P(\text{Tu}, -1 \mid -1), \ldots, P(\text{Fr}, 1 \mid 1).\]

Если, допустим, известны все вероятности, стоящие в формуле Байеса справа от знака равенства, тогда можно рассчитать условную вероятность слева. Тогда по любому набору значений признаков можно найти вероятность отнести предъявленный объект к каждому из классов. Какая вероятность окажется больше, тот класс мы и выберем. Вероятность в знаменателе формулы Байеса можно не рассматривать, поскольку она не зависит от того, к какому классу принадлежит тот или иной объект. Чтобы выбрать наибольшее значение из всех рассчитанных, нам достаточно знать только вероятности, стоящие в числителе.

Байесовский классификатор требует обучения: мы предъявляем ему набор объектов с различными значениями признаков (в нашем примере это изменения цен за какой-либо отрезок истории) и сообщаем, к какому классу должен быть отнесён каждый из этих объектов. Описанный процесс называется обучением с учителем (учитель – это тот, кто знает, к какому классу должны быть отнесены объекты). В результате обучения классификатор должен найти значения всех вероятностей, стоящих в формуле Байеса справа от знака равенства в числителе.

Но задачу можно упростить, если предположить, что все признаки независимы друг от друга. Тогда можно вероятность сложного события записать в виде произведения вероятностей каждого события в отдельности:

\[\mathbb{P}(X_1, X_2,\ldots, X_n \mid C) = \mathbb{P}(X_1 \mid C) \cdot \mathbb{P}(X_2 \mid C)\cdot\ldots\cdot \mathbb{P}(X_n \mid C).\]

Что это даёт? Намного проще найти \(n\) одномерных функций, чем одну \(n\)-мерную функцию. Для нашего примера вместо 45-ти условных вероятностей \(\mathbb{P}(X_1, X_2 \mid C)\) надо оценить по экспериментальным данным из обучающей выборки независимо друг от друга 15 значений \(\mathbb{P}(X_1 \mid C)\) (5 рабочих дней недели на 3 состояния рынка) и 9 значений \(\mathbb{P}(X_2 \mid C)\) (три состояния рынка вчера на три состояния сегодня), т.е. всего 24 искомых параметров вместо 45. При увеличении размерности задачи (количества признаков) выигрыш ещё более заметен.

Таким образом, наи́вный ба́йесовский классифика́тор — это простой вероятностный классификатор, основанный на применении Теоремы Байеса с (наивными) предположениями о независимости предикторов:

\[\mathbb{P}(C \mid X_1, X_2,\ldots, X_n) = \frac{\mathbb{P}(X_1 \mid C) \cdot \mathbb{P}(X_2 \mid C)\cdot\ldots\cdot \mathbb{P}(X_n \mid C) \cdot\mathbb{P}(C) }{\mathbb{P}(X_1)\cdot \mathbb{P}(X_2)\cdot\ldots\cdot \mathbb{P}(X_n)}.\]

Реализация наивного байесовского классификатора на R

Теперь рассмотрим работу наивного байесовского классификатора на реальных биржевых данных.

Создадим новый R-скрипт или будем выполнять следующие команды в командной строке R:

library("quantmod")  # Для скачивания котировок.
library("lubridate") # Для работы с датами.
library("e1071")     # Алгоритмы машинного обучения.
startDate = as.Date("2009-01-01") # Начальная дата.
endDate = as.Date("2013-12-31")   # Конечная дата.

# Скачиваем котировки индекса S&P500 с Yahoo Finance:
stock = getSymbols("^GSPC", src = "yahoo", from = startDate, to = endDate, auto.assign=F)
colnames(stock) = c("Open", "High", "Low", "Close", "Volume", "Adjusted") # Имена столбцов.

dayOfWeek = wday(stock, label=T, abbr=T) # День недели для каждой даты.
dr = dailyReturn(Ad(stock), type='arithmetic')  # Относительные приращения цены за каждый день.
# Удалим первый элемент в обоих наборах данных (для него нет предыдущего дня):
dr = dr[-1];  dayOfWeek = dayOfWeek[-1]

# Максимальное, минимальное и среднее значения относительного приращения цен:
max(dr); min(dr); mean(dr)

Получим диапазон колебаний относительного приращения цен закрытия (simple returns) и среднее значение:

[1] 0.07074453
[1] -0.06663443
[1] 0.0006184646

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

s = 0.0035 # Максимальное относительное приращение цены, которое отличает флет от тренда.
# Если относительное приращение цены больше s, то имеем тренд, иначе флет:
classes = ifelse(dr > s, "1", ifelse(dr < -s, "-1", "0"))

# Набор данных с двумя предикторами и столбцом значений классов:
dataset2 = data.frame(classes[-1], dayOfWeek[-1], classes[-length(classes)])
colnames(dataset2) = c("Класс", "ДеньНедели", "НаправлениеВчера")

Разделим наш набор данных на две неравные части. Три четверти всех дней будем использовать для обучения классификатора, а последнюю четверть — для проверки его работы:

n = nrow(dataset2) # Длина всего набора данных (кол-во дней)
ntrain = round(3 * n / 4, digits=0) # Первые три четверти - обучающие.
train2 = dataset2[1:ntrain-1,] # Получили обучающую выборку из исходного набора данных.
test2 = dataset2[ntrain:n,]  # Получили проверочную выборку из исходного набора данных.

Теперь построим наивный байесовский классификатор и передадим ему данные для обучения:

# Обучающий набор предикторов - столбцы 2 и 3, значения классов - столбец 1:
model2 = naiveBayes(train2[,2:3], train2[,1])
model2 # Вывели параметры модели на экран.

Получили априорные вероятности состояний рынка \(\mathbb{P}(C)\):

A-priori probabilities:
train2[, 1]
        0          1         -1
0.3319194  0.3753977  0.2926829

а также апостериорные условные вероятности \(\mathbb{P}(X_1 \mid C)\):

Conditional probabilities:
                ДеньНедели
train2[, 1]     Sun        Mon       Tues        Wed      Thurs        Fri        Sat
   0      0.0000000  0.1980831  0.2076677  0.2108626  0.1853035  0.1980831  0.0000000
   1      0.0000000  0.1723164  0.2062147  0.2033898  0.2288136  0.1892655  0.0000000
  -1      0.0000000  0.1884058  0.2065217  0.2028986  0.1920290  0.2101449  0.0000000

и \(\mathbb{P}(X_2 \mid C)\):

             НаправлениеВчера
train2[, 1]       0          1         -1
   0      0.3514377  0.4249201  0.2236422
   1      0.3079096  0.3502825  0.3418079
  -1      0.3405797  0.3514493  0.3079710

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

pr2train = predict(model2, train2[,2:3]) # Предсказанные состояния рынка для обучающей выборки.
pr2test = predict(model2, test2[,2:3]) # Предсказанные состояния рынка для проверочной выборки.

# Таблица правильности предсказаний для обучающей выборки:
t2train = table(pr2train, train2[,1], dnn=list('Предсказано', 'На самом деле'))
t2train # Вывели таблицу на экран

# Таблица правильности предсказаний для проверочной выборки:
t2test = table(pr2test, test2[,1], dnn=list('Предсказано', 'На самом деле'))
t2test # Вывели таблицу на экран

Для обучающей выборки (943 дня) получим:

            На самом деле
Предсказано   0    1   -1
        0   201  175  153
        1   112  179  123
       -1     0    0    0

Это значит, что из всех случаев, когда классификатор предсказал тренд вверх (таких дней было 112+179+123=414), в 112-ти случаях на самом деле наблюдался флет, в 179-ти случаях действительно был тренд вверх и в 123-х случаях рынок падал. Нисходящий тренд не был предсказан ни разу.

Для проверочной выборки (316 дней), т.е. по тем данным, которые не использовались при обучении:

         На самом деле
Предсказано  0   1  -1
        0   85  65  51
        1   41  50  24
       -1    0   0   0

Таким образом, восходящий тренд был предсказан 41+50+24=115 раз (в среднем каждый третий день), из них в 41-м случае на самом деле наблюдался флет, в 50-ти случаях мы не ошиблись, а в 24-х случаях тренд оказался нисходящим. Нисходящий тренд не был предсказан ни разу.

Допустим, в те дни, когда предсказан флет, мы не открываем позиций по данному инструменту; а в дни, когда предсказан восходящий тренд, — открываем длинную позицию. Если мы открыли позицию, а на самом деле день оказался флетовым, то мы, скорее всего, ничего не заработали, но и ничего не потеряли. Если направление угадано правильно, мы в конце дня закрываем позицию с прибылью. Если мы ошиблись с направлением, то фиксируем убыток. Получается, что из 115-ти дней, когда мы открывали длинные позиции, согласно прогнозу классификатора, 50 дней оказались прибыльными, а 24 дня — убыточными. Не такой уж плохой результат: 50/24 = 2,08.

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

Пересечение скользящих средних

. . . . . . .

Список литературы

. . . . . . .

Ссылки



Комментарии

Комментариев пока нет.

* Обязательные поля
(Не публикуется)
 
Жирный Курсив Подчеркнутый Перечеркнутый Степень Индекс Код PHP Код Кавычки Вставить линию Вставить маркированный список Вставить нумерованный список Вставить ссылку Вставить e-mail Вставить изображение Вставить видео
 
Улыбка Печаль Удивление Смех Злость Язык Возмущение Ухмылка Подмигнуть Испуг Круто Скука Смущение Несерьёзно Шокирован
 
1000
Captcha
Refresh
 
Введите код:
 
Запомнить информацию введенную в поля формы.