Полёты на падающих снежинках02.12.2017 18:33

Усложним задачу. Пусть количество предикторов будет произвольным. Это уже не прямая, а гиперплоскость в многомерном пространстве. И сменим тип задачи с аппроксимации на классификацию. Например, с помощью логистической регрессии была решена задача бинарной классификации. Мы хотим использовать эту статистическую модель в своём невероятно сильно нагруженном сервисе, скажем, для классификации пользователей. Для решения задачи нам не нужен сам алгоритм обучения модели, а только параметры разделяющей гиперплоскости.

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

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

# Загрузка данных
dataset <- read.csv('dataset.csv')

# Для наглядности визуально отобразим данные
pairs(dataset, col = factor(dataset$class))

# Тренировка модели логистической регрессии
model <- glm(formula = class ~ ., data = dataset, family = binomial)

# Получение необходимых коэффициентов
b <- model$coefficients

# Вероятности нам не нужны, просто разделеним точки
nc <- (b[1] + (dataset$alpha * b[2]) + (dataset$beta * b[3])) > 0
plot(dataset$alpha, dataset$beta, col = factor(nc))

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

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

test <- factor(as.logical(dataset$class))

length(nc)
# 345

table(test == nc)
# FALSE  TRUE
#    17   328

328 / 345
# 0.9507246

А что, если доля наблюдений одного из классов будет всего одна тысячная процента? Тут даже выдающий константу классификатор покажет замечательный результат. Есть смысл отобразить confusion matrix. В бинарной классификации могут быть только четыре исхода предсказания метки класса. Истинным положительным результатом называют правильное угадывание положительного класса (TRUE предсказан как TRUE). Истинным отрицательным — верное угадывание отрицательного класса (FALSE предсказан как FALSE). Логично предположить, что ложный положительный — это FALSE предсказанный как TRUE, а ложный отрицательный — это TRUE предсказанный как FALSE. Посмотрим на конкретном примере:

table(nc, test)
#          test
#   nc   FALSE TRUE
#  FALSE   120    8
#  TRUE      9  208

Где «nc» — это ответ классификатора, а «test» — истинный ответ. Определим долю положительных ответов классификатора, которые действительно были положительными (точность, precision). А также полноту (recall), т.е. какую долю настоящих положительных удалось выявить этой моделью. Смысл у этих показателей следующий: точность показывает степень доверия классификатору, если он дал положительный ответ. Другими словами, насколько мы можем быть уверены, что это действительно положительный класс. А вот полнота показывает охват выявляющей способности, т.е. долю выявленных положительных. Если мы боимся ошибочно назвать класс положительным, то важнее точность. Когда нужно найти как можно больше положительных, то важнее полнота.

library(caret)

precision <- posPredValue(factor(nc), test, positive = T)
# 0.9585253

recall <- sensitivity(factor(nc), test, positive = T)
# 0.962963

# Метрика F1 является средним гармоническим полноты и точности
f1 <- (2 * precision * recall) / (precision + recall)
# 0.960739

Рассчитаем вручную точность и полноту. Достаточно подставить значения из confusion matrix:

# precision (ИстП / ИстП + ЛожП)
208 / (208 + 9)
# 0.9585253

# recall (ИстП / ИстП + ЛожОтр)
208 / (208 + 8)
# 0.962963

4. Оценка важности предикторов

Допустим, производилось комплексное психологическое исследование людей. Одна половина испытуемых страдает неврозом, а у другой всё хорошо. В чём отличия? Или иной пример: измерили показатели оборудования — одни устройства работают хорошо, а с другими проблема. Что влияет на это? Интуитивно можно предположить, что следует попробовать поискать корреляцию между зависимой переменной и каждым из предикторов. Вдруг окажется, что, скажем, метрика качества топлива сильно коррелирует с долговечностью (чем лучше топливо, тем долговечнее работа устройства). Или посмотреть различия средних значений предикторов у наблюдений разных классов.

Однако, попробуем визуально отобразить наш набор данных. Видна некоторая условная граница, после которой точки меняют цвет. Она проходит в районе 0.5 по предиктору «alpha». Это видно и по гистограмме его распределения (разными цветами показано). А по предиктору «beta» таких явных различий не наблюдается.

Интуитивно понятно, что несмотря на ошибку (некоторые точки окажутся неправильно классифицированными) такой критерий будет эффективнее решать проблему, чем любое возможное разделение по предиктору «beta». Следовательно, важность «alpha» будет достаточно высокой. После разделения вероятность встретить точку нужного класса значительно возрастает. Чтобы выявить разницу вначале вычислим вероятность без разделения.

table(dataset$class)
#   0   1 
# 129 216

length(dataset$class)
# 345

c(129 / 345, 216 / 345)
# 0.373913 0.626087

Зная вероятности вычислим метрику однородности данных (если только представители одного класса, то показатель Gini impurity будет равен 0). Gini impurity вычисляется как сумма квадратов вероятностей, которую отняли от единицы.

# Gini impurity
gini <- function(p) {
    (1 - sum(p ^ 2)) 
}

gini(c(0.373913, 0.626087))
# 0.4682041

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

node_1 <- subset(dataset, alpha > .5)

table(node_1$class)
#   0   1 
#  25 197

length(node_1$class)
# 222

# Разнообразие снижено
gini(c(25/222, 197/222))
# 0.199862

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

node_2 <- subset(node_1, beta < .29)

table(node_2$class)
#   1
#  34

length(node_2$class)
# 34

# Это уже будет терминальный узел
gini(c(0/34, 34/34))
# 0

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

На этом наборе данных было максимально очевидно, что один из предикторов обладает очень высокой значимостью. Посмотрим, как он коррелирует с меткой класса:

Видны различия и в среднем значении разных классов. Предлагаю посмотреть описательную статистику (примеры на Python и R):

Оценка важности предикторов посредством Random Forest:

5. Выводы

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

6. Приложение

Использованные фрагменты исходного кода (Python) для выявления степени важности предикторов:

import pandas as pd

dataset = pd.read_csv('dataset.csv')
dataset.info()
dataset.sample(5)
dataset.describe()
dataset.corr()
dataset.groupby('class').mean()

import seaborn as sns
import matplotlib.pyplot as plt

dataset.hist(bins = 20, figsize = (6, 6))
plt.show()

sns.heatmap(dataset.corr(), square = True, annot = True)
plt.show()

sns.pairplot(
    data = dataset,
    hue = 'class',
    size = 2,
    palette = 'seismic'
)
plt.show()

import pandas as pd
from sklearn.metrics import classification_report
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier

dataset = pd.read_csv('dataset.csv')

X_train, X_test, y_train, y_test = train_test_split(
    dataset,
    dataset.pop('class'),
    test_size = .5
)

model = RandomForestClassifier(
    n_estimators = 1000,
    max_depth = 100, 
    max_features = 2
).fit(X_train, y_train)

print(classification_report(y_test, model.predict(X_test)))

pd.Series(model.feature_importances_, index = X_train.columns)

from catboost import CatBoostClassifier
from sklearn.model_selection import KFold
from sklearn.metrics import confusion_matrix, f1_score

kf = KFold(n_splits = 3, shuffle = True)
dataset = pd.read_csv('dataset.csv')

y = dataset.pop('class').values
X = dataset.values
columns = dataset.columns

for train_index, test_index in kf.split(X):
    X_train, X_test = X[train_index], X[test_index]
    y_train, y_test = y[train_index], y[test_index]

    model = CatBoostClassifier(
        calc_feature_importance = True
    ).fit(X_train, y_train)

    y_pred = model.predict(X_test)

    print(confusion_matrix(y_test, y_pred))
    print(f1_score(y_test, y_pred))
    print(list(zip(columns, model.feature_importances_)))

© Habrahabr.ru