Хочу поделиться опытом работы с задачей известного конкурса по машинному обучению от Kaggle. Этот конкурс позиционируется как конкурс для начинающих, а у меня как раз не было почти никакого практического опыта в этой области. Я немного знал теорию, но с реальными данными дела почти не имел и с питоном плотно не работал. В итоге, потратив пару предновогодних вечеров, набрал 0.80383 (первая четверть рейтинга).



В общем эта статья для еще начинающих от уже начавшего.



Титаник



Включаем подходящую для работы музыку и начинаем исследование.



Конкурс «по титанику» был уже неоднократно отмечен на Хабре. Особенно хотелось бы отметить последнюю статью из списка — оказывается, что исследование данных может быть не менее интригующим чем хороший детективный роман, а отличный результат (0.81340) можно получить не только Random Forest классификатором.



Также хотелось бы отметить статью о другом конкурсе. Из нее можно понять каким именно образом должен работать мозг исследователя и что большая часть времени должна быть уделена предварительному анализу и обработке данных.



Инструментарий



Для решения задачи я использую Python-стек технологий. Этот подход не является единственно возможным: существуют R, Matlab, Mathematica, Azure Machine Learning, Apache Weka, Java-ML и я думаю список можно еще долго продолжать. Использование Python имеет ряд преимуществ: библиотек действительно много и они отличного качества, а поскольку большинство из них представляют собой обертки над C-кодом они еще и достаточно быстрые. К тому же построенная модель может быть легко запущена в эксплуатацию.

Должен признаться, что я не очень большой поклонник скриптовых нестрого-типизированных языков, но богатство библиотек для python не позволяет его как-либо игнорировать.

Запускать все будем под Linux (Ubuntu 14.04). Понадобятся: python 2.7, seaborn, matplotlib, sklearn, xgboost, pandas. Вообще обязательны только pandas и sklearn, а остальные нужны для иллюстрации.

Под Linux библиотеки для Python можно устанавливать двумя способами: штатным пакетным (deb) менеджером или через питоновскую утилиту pip.

Установка deb-пакетами проще и быстрее, но зачастую библиотеки там устаревшие (стабильность превыше всего).

# Установка будет произведена в /usr/lib/python2.7/dist-packages/
$ sudo apt-get install python-matplotlib


Установка пакетов через pip дольше (будет произведена компиляция), но с ней можно рассчитывать на получение свежих версий пакетов.

# Установка будет произведена в /usr/local/lib/python2.7/dist-packages/
$ sudo pip install matplotlib


Так каким именно образом лучше устанавливать пакеты? Я использую компромисс: массивные и требующие множества зависимостей для сборки NumPy и SciPy ставлю из DEB-пакетов.

$ sudo apt-get install python 
$ sudo apt-get install python-pip
$ sudo apt-get install python-numpy
$ sudo apt-get install python-scipy
$ sudo apt-get install ipython


А остальные, более легкие пакеты, устанавливаю через pip.

$ sudo pip install pandas 
$ sudo pip install matplotlib==1.4.3
$ sudo pip install skimage
$ sudo pip install sklearn
$ sudo pip install seaborn
$ sudo pip install statsmodels
$ sudo pip install xgboost


Если я что-то забыл — то все нужные пакеты как правило легко вычисляются и устанавливаются аналогичным способом.

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



Данные



Скачаем исходные данные для задачи и посмотрим что же нам выдали.

$ wc -l train.csv test.csv 
  892 train.csv
  419 test.csv
 1311 total


Данных прямо скажем не очень много — всего 891 пассажир в train-выборке, и 418 в test-выборке (одна строчка идет на заголовок со списком полей).

Откроем train.csv в любом табличном процессоре (я использую LibreOffice Calc) чтобы визуально посмотреть данные.

$ libreoffice --calc train.csv


Видим следующее:
  • Возраст заполнен не у всех
  • Билеты имеют какой-то странный и неконсистентный формат
  • В именах есть титул (мисс, мистер, миссис и т.д.)
  • Номера кают прописаны мало у кого (тут есть леденящая душу история почему)
  • В существующих номерах кают видимо прописан код палубы (так и есть, как оказалось)
  • Также, согласно статье, в номере каюты зашифрована сторона
  • Сортируем по имени. Видно, что многие путешествовали семьями, и виден масштаб трагедии — часто семьи оказались разделены, выжила лишь часть
  • Сортируем по билету. Видно, что по одному коду билета путешествовало сразу несколько людей, причем часто — с разными фамилиями. Беглый взгляд вроде бы показывает, что люди с одинаковым номером билета часто разделяют одну и ту же участь.
  • У части пассажиров не проставлен порт посадки


Вроде приблизительно все ясно, переходим непостредственно к работе с данными.

Загрузка данных



Чтобы не зашумлять дальшейший код сразу приведу все используемые импорты:
Шапка скрипта
# coding=utf8

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import xgboost as xgb
import re
import seaborn as sns
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import LogisticRegression
from sklearn.linear_model import SGDClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.grid_search import GridSearchCV
from sklearn.neighbors import KNeighborsClassifier
from sklearn.naive_bayes import GaussianNB
from sklearn.svm import SVC
from sklearn.feature_selection import SelectKBest, f_classif
from sklearn.cross_validation import StratifiedKFold
from sklearn.cross_validation import KFold
from sklearn.cross_validation import cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn import metrics

pd.set_option('display.width', 256)


Наверное большая часть моей мотивации к написанию этой статьи вызвана восторгом от работы с пакетом pandas. Я знал о существовании этой технологии, но не мог даже представить насколько приятно с ней работать. Pandas — это Excel в командной строке с удобной функциональностью по вводу/выводу и обработке табличных данных.

Загружаем обе выборки.

train_data = pd.read_csv("data/train.csv")
test_data = pd.read_csv("data/test.csv")


Собираем обе выборки (train-выборку и test-выборку) в одну суммарную all-выборку.

all_data = pd.concat([train_data, test_data])


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

Анализ данных



Анализ данных в Python можно осуществлять сразу несколькими способами, например:

  • Вручную готовить данные и выводить через matplotlib
  • Использовать все готовое в seaborn
  • Использовать текстовый вывод с группировкой в pandas


Мы попробуем все три, но сначала запустим самый простой текстовый вариант. Выведем статистику выживаемости в зависимости от класса и пола.

print("===== survived by class and sex")
print(train_data.groupby(["Pclass", "Sex"])["Survived"].value_counts(normalize=True))

Результат
===== survived by class and sex
Pclass  Sex     Survived
1       female  1           0.968085
                0           0.031915
        male    0           0.631148
                1           0.368852
2       female  1           0.921053
                0           0.078947
        male    0           0.842593
                1           0.157407
3       female  0           0.500000
                1           0.500000
        male    0           0.864553
                1           0.135447
dtype: float64



Видим, что в шлюпки действительно сажали сначала женщин — шанс женщины на выживаемость составляет 96.8%, 92.1% и 50% в зависимости от класса билета. Шанс на выживание мужчины гораздно ниже и составляет соответственно 36.9%, 15.7% и 13.5%.

С помощью pandas быстро посчитаем сводку по всем числовым полям обеих выборок — отдельно по мужчинам и по женщинам.

describe_fields = ["Age", "Fare", "Pclass", "SibSp", "Parch"]

print("===== train: males")
print(train_data[train_data["Sex"] == "male"][describe_fields].describe())

print("===== test: males")
print(test_data[test_data["Sex"] == "male"][describe_fields].describe())

print("===== train: females")
print(train_data[train_data["Sex"] == "female"][describe_fields].describe())

print("===== test: females")
print(test_data[test_data["Sex"] == "female"][describe_fields].describe())

Результат
===== train: males
              Age        Fare      Pclass       SibSp       Parch
count  453.000000  577.000000  577.000000  577.000000  577.000000
mean    30.726645   25.523893    2.389948    0.429809    0.235702
std     14.678201   43.138263    0.813580    1.061811    0.612294
min      0.420000    0.000000    1.000000    0.000000    0.000000
25%     21.000000    7.895800    2.000000    0.000000    0.000000
50%     29.000000   10.500000    3.000000    0.000000    0.000000
75%     39.000000   26.550000    3.000000    0.000000    0.000000
max     80.000000  512.329200    3.000000    8.000000    5.000000
===== test: males
              Age        Fare      Pclass       SibSp       Parch
count  205.000000  265.000000  266.000000  266.000000  266.000000
mean    30.272732   27.527877    2.334586    0.379699    0.274436
std     13.389528   41.079423    0.808497    0.843735    0.883745
min      0.330000    0.000000    1.000000    0.000000    0.000000
25%     22.000000    7.854200    2.000000    0.000000    0.000000
50%     27.000000   13.000000    3.000000    0.000000    0.000000
75%     40.000000   26.550000    3.000000    1.000000    0.000000
max     67.000000  262.375000    3.000000    8.000000    9.000000
===== train: females
              Age        Fare      Pclass       SibSp       Parch
count  261.000000  314.000000  314.000000  314.000000  314.000000
mean    27.915709   44.479818    2.159236    0.694268    0.649682
std     14.110146   57.997698    0.857290    1.156520    1.022846
min      0.750000    6.750000    1.000000    0.000000    0.000000
25%     18.000000   12.071875    1.000000    0.000000    0.000000
50%     27.000000   23.000000    2.000000    0.000000    0.000000
75%     37.000000   55.000000    3.000000    1.000000    1.000000
max     63.000000  512.329200    3.000000    8.000000    6.000000
===== test: females
              Age        Fare      Pclass       SibSp       Parch
count  127.000000  152.000000  152.000000  152.000000  152.000000
mean    30.272362   49.747699    2.144737    0.565789    0.598684
std     15.428613   73.108716    0.887051    0.974313    1.105434
min      0.170000    6.950000    1.000000    0.000000    0.000000
25%     20.500000    8.626050    1.000000    0.000000    0.000000
50%     27.000000   21.512500    2.000000    0.000000    0.000000
75%     38.500000   55.441700    3.000000    1.000000    1.000000
max     76.000000  512.329200    3.000000    8.000000    9.000000



Видно, что по средним и перцентилям все в целом ровненько. Но у мужчин по выборкам отличаются максимумы по возрасту и по стоимости билета. У женщин в обеих выборках также есть различие в максимуме возраста.

Сборка дайджеста по данным



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

class DataDigest:

    def __init__(self):
        self.ages = None
        self.fares = None
        self.titles = None
        self.cabins = None
        self.families = None
        self.tickets = None

def get_title(name):
    if pd.isnull(name):
        return "Null"

    title_search = re.search(' ([A-Za-z]+)\.', name)
    if title_search:
        return title_search.group(1).lower()
    else:
        return "None"


def get_family(row):
    last_name = row["Name"].split(",")[0]
    if last_name:
        family_size = 1 + row["Parch"] + row["SibSp"]
        if family_size > 3:
            return "{0}_{1}".format(last_name.lower(), family_size)
        else:
            return "nofamily"
    else:
        return "unknown"


data_digest = DataDigest()
data_digest.ages = all_data.groupby("Sex")["Age"].median()
data_digest.fares = all_data.groupby("Pclass")["Fare"].median()
data_digest.titles = pd.Index(test_data["Name"].apply(get_title).unique())
data_digest.families = pd.Index(test_data.apply(get_family, axis=1).unique())
data_digest.cabins = pd.Index(test_data["Cabin"].fillna("unknown").unique())
data_digest.tickets = pd.Index(test_data["Ticket"].fillna("unknown").unique())


Небольшое пояснение по полям дайджеста:
  • ages — справочник медиан возрастов в зависимости от пола;
  • fares — справочник медиан стоимости билетов в зависимости от класса билета;
  • titles — справочник титулов;
  • families — справочник идентификаторов семейств (фамилия + количество членов семьи);
  • cabins — справочник идентификаторов кают;
  • tickets — справочник идентификаторов билетов.


Справочники для восстановления отсутствующих данных (медианы) мы строим по комбинированной выборке. А вот справочники для перевода категориальных признаков — только по тестовым данным. Идея заключалась в следующем: допустим в train-наборе у нас есть фамилия «Иванов», а в test-наборе этой фамилии нет. Знание внутри классификатора о том что «Иванов» выжил (или не выжил) никак не поможет в оценке test-набора, поскольку в test-наборе этой фамилии все равно нет. Поэтому в справочник добавляем только те фамилии которые есть в test-наборе. Еще более правильным способом будет добавлять в справочник только переcечение признаков (только те признаки которые есть в обоих наборах) — я попробовал, но результат верификации ухудшился на 3 процента.

Выделяем признаки



Теперь нам нужно выделить признаки. Как уже было сказано — многие классификаторы умеют работать только с числами, поэтому нам нужно:

  • Перевести категории в числовое представление
  • Выделить неявные признаки, то есть те которые явно не заданы (титул, палуба)
  • Что-то сделать с отсутствующими значениями


Есть два способа преобразования категориального признака в числовой. Можем рассмотреть задачу на примере пола пассажира.

В первом варианте мы просто меняем пол на некоторое число, например мы можем заменить female на 0, а male на 1 (кругляшок и палочка — очень удобно запоминать). Такой вариант не увеличивает число признаков, однако внутри признака для его значений теперь появляется отношение «больше» и «меньше». В случае когда значений много, такое неожиданное свойство признака не всегда желательно и может привести к проблемам в геометрических классификаторах.

Второй вариант преобразования — завести две колонки «sex_male» и «sex_female». В случае мужского пола мы будем присваивать sex_male=1, sex_female=0. В случае женского пола наоборот: sex_male=0, sex_female=1. Отношений «больше»/«меньше» мы теперь избегаем, однако теперь у нас стало больше признаков, а чем больше признаков тем больше данных для тренировки классификатора нам нужно — эта проблема известна как «проклятие размерности». Особенно сложной становится ситуация когда значений признаков очень много, например идентификаторы билетов, в таких случаях можно например откинуть редко встречающиеся значения подставив вместо них какой-нибудь специальный тег — таким образом уменьшив итоговое количество признаков после расширения.

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

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

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

Создаем метод для преобразования наборов данных.

def get_index(item, index):
    if pd.isnull(item):
        return -1

    try:
        return index.get_loc(item)
    except KeyError:
        return -1


def munge_data(data, digest):
    # Age - замена пропусков на медиану в зависимости от пола
    data["AgeF"] = data.apply(lambda r: digest.ages[r["Sex"]] if pd.isnull(r["Age"]) else r["Age"], axis=1)

    # Fare - замена пропусков на медиану в зависимости от класса
    data["FareF"] = data.apply(lambda r: digest.fares[r["Pclass"]] if pd.isnull(r["Fare"]) else r["Fare"], axis=1)

    # Gender - замена
    genders = {"male": 1, "female": 0}
    data["SexF"] = data["Sex"].apply(lambda s: genders.get(s))

    # Gender - расширение
    gender_dummies = pd.get_dummies(data["Sex"], prefix="SexD", dummy_na=False)
    data = pd.concat([data, gender_dummies], axis=1)

    # Embarkment - замена
    embarkments = {"U": 0, "S": 1, "C": 2, "Q": 3}
    data["EmbarkedF"] = data["Embarked"].fillna("U").apply(lambda e: embarkments.get(e))

    # Embarkment - расширение
    embarkment_dummies = pd.get_dummies(data["Embarked"], prefix="EmbarkedD", dummy_na=False)
    data = pd.concat([data, embarkment_dummies], axis=1)

    # Количество родственников на борту
    data["RelativesF"] = data["Parch"] + data["SibSp"]

    # Человек-одиночка?
    data["SingleF"] = data["RelativesF"].apply(lambda r: 1 if r == 0 else 0)

    # Deck - замена
    decks = {"U": 0, "A": 1, "B": 2, "C": 3, "D": 4, "E": 5, "F": 6, "G": 7, "T": 8}
    data["DeckF"] = data["Cabin"].fillna("U").apply(lambda c: decks.get(c[0], -1))

    # Deck - расширение
    deck_dummies = pd.get_dummies(data["Cabin"].fillna("U").apply(lambda c: c[0]), prefix="DeckD", dummy_na=False)
    data = pd.concat([data, deck_dummies], axis=1)

    # Titles - расширение
    title_dummies = pd.get_dummies(data["Name"].apply(lambda n: get_title(n)), prefix="TitleD", dummy_na=False)
    data = pd.concat([data, title_dummies], axis=1)

    # амена текстов на индекс из соответствующего справочника или -1 если значения в справочнике нет (расширять не будем)
    data["CabinF"] = data["Cabin"].fillna("unknown").apply(lambda c: get_index(c, digest.cabins))

    data["TitleF"] = data["Name"].apply(lambda n: get_index(get_title(n), digest.titles))

    data["TicketF"] = data["Ticket"].apply(lambda t: get_index(t, digest.tickets))

    data["FamilyF"] = data.apply(lambda r: get_index(get_family(r), digest.families), axis=1)

    # для статистики
    age_bins = [0, 5, 10, 15, 20, 25, 30, 40, 50, 60, 70, 80, 90]
    data["AgeR"] = pd.cut(data["Age"].fillna(-1), bins=age_bins).astype(object)

    return data


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


В общем, мы добавляем в признаки вообще все что приходит в голову. Видно что некоторые признаки дублируют друг друга (например расширение и замена пола), некоторые явно коррелируют друг с другом (класс билета и стоимость билета), некоторые явно бессмысленны (вряд ли порт посадки влияет на выживаемость). Разбираться со всем этим будем позже — когда будем производить отбор признаков для обучения.

Преобразуем оба имеющихся набора и также опять создаем объединенный набор.

train_data_munged = munge_data(train_data, data_digest)
test_data_munged = munge_data(test_data, data_digest)
all_data_munged = pd.concat([train_data_munged, test_data_munged])


Хотя мы нацелились на использование Random Forest хочется попробовать и другие классификаторы. А с ними есть следующеая проблема: многие классификаторы чувствительны к масштабу признаков. Другими словами если у нас есть один признак со значениями от [-10,5] и второй признак со значениями [0,10000] то одинаковая в процентном отношении ошибка на обоих признаках будет приводить к большому отличию в абсолютном значении и классификатор будет трактовать второй признак как более важный.

Чтобы избежать этого мы приводим все числовые (а других у нас уже нет) признаки к одинаковой шкале [-1,1] и нулевому среднему значению. Сделать это в sklearn можно очень просто.

scaler = StandardScaler()
scaler.fit(all_data_munged[predictors])

train_data_scaled = scaler.transform(train_data_munged[predictors])
test_data_scaled = scaler.transform(test_data_munged[predictors])


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

Выбор признаков



Ну и пришел момент когда мы можем отобрать те признаки, с которыми будем работать дальше.

predictors = ["Pclass",
              "AgeF",
              "TitleF",
              "TitleD_mr", "TitleD_mrs", "TitleD_miss", "TitleD_master", "TitleD_ms", 
              "TitleD_col", "TitleD_rev", "TitleD_dr",
              "CabinF",
              "DeckF",
              "DeckD_U", "DeckD_A", "DeckD_B", "DeckD_C", "DeckD_D", "DeckD_E", "DeckD_F", "DeckD_G",
              "FamilyF",
              "TicketF",
              "SexF",
              "SexD_male", "SexD_female",
              "EmbarkedF",
              "EmbarkedD_S", "EmbarkedD_C", "EmbarkedD_Q",
              "FareF",
              "SibSp", "Parch",
              "RelativesF",
              "SingleF"]


Просто ставим комментарий на ненужных и запускаем обучение. Какие именно не нужны — решать вам.

Еще раз анализ



Поскольку теперь у нас появилась колонка в которой прописан диапазон в который попадает возраст пассажира — оценим выживаемость в зависимости от возраста (диапазона).

print("===== survived by age")
print(train_data.groupby(["AgeR"])["Survived"].value_counts(normalize=True))

print("===== survived by gender and age")
print(train_data.groupby(["Sex", "AgeR"])["Survived"].value_counts(normalize=True))

print("===== survived by class and age")
print(train_data.groupby(["Pclass", "AgeR"])["Survived"].value_counts(normalize=True))

Результат
===== survived by age
AgeR      Survived
(0, 5]    1           0.704545
          0           0.295455
(10, 15]  1           0.578947
          0           0.421053
(15, 20]  0           0.656250
          1           0.343750
(20, 25]  0           0.655738
          1           0.344262
(25, 30]  0           0.611111
          1           0.388889
(30, 40]  0           0.554839
          1           0.445161
(40, 50]  0           0.616279
          1           0.383721
(5, 10]   0           0.650000
          1           0.350000
(50, 60]  0           0.595238
          1           0.404762
(60, 70]  0           0.764706
          1           0.235294
(70, 80]  0           0.800000
          1           0.200000
dtype: float64
===== survived by gender and age
Sex     AgeR      Survived
female  (0, 5]    1           0.761905
                  0           0.238095
        (10, 15]  1           0.750000
                  0           0.250000
        (15, 20]  1           0.735294
                  0           0.264706
        (20, 25]  1           0.755556
                  0           0.244444
        (25, 30]  1           0.750000
                  0           0.250000
        (30, 40]  1           0.836364
                  0           0.163636
        (40, 50]  1           0.677419
                  0           0.322581
        (5, 10]   0           0.700000
                  1           0.300000
        (50, 60]  1           0.928571
                  0           0.071429
        (60, 70]  1           1.000000
male    (0, 5]    1           0.652174
                  0           0.347826
        (10, 15]  0           0.714286
                  1           0.285714
        (15, 20]  0           0.870968
                  1           0.129032
        (20, 25]  0           0.896104
                  1           0.103896
        (25, 30]  0           0.791667
                  1           0.208333
        (30, 40]  0           0.770000
                  1           0.230000
        (40, 50]  0           0.781818
                  1           0.218182
        (5, 10]   0           0.600000
                  1           0.400000
        (50, 60]  0           0.857143
                  1           0.142857
        (60, 70]  0           0.928571
                  1           0.071429
        (70, 80]  0           0.800000
                  1           0.200000
dtype: float64
===== survived by class and age
Pclass  AgeR      Survived
1       (0, 5]    1           0.666667
                  0           0.333333
        (10, 15]  1           1.000000
        (15, 20]  1           0.800000
                  0           0.200000
        (20, 25]  1           0.761905
                  0           0.238095
        (25, 30]  1           0.684211
                  0           0.315789
        (30, 40]  1           0.755102
                  0           0.244898
        (40, 50]  1           0.567568
                  0           0.432432
        (50, 60]  1           0.600000
                  0           0.400000
        (60, 70]  0           0.818182
                  1           0.181818
        (70, 80]  0           0.666667
                  1           0.333333
2       (0, 5]    1           1.000000
        (10, 15]  1           1.000000
        (15, 20]  0           0.562500
                  1           0.437500
        (20, 25]  0           0.600000
                  1           0.400000
        (25, 30]  0           0.580645
                  1           0.419355
        (30, 40]  0           0.558140
                  1           0.441860
        (40, 50]  1           0.526316
                  0           0.473684
        (5, 10]   1           1.000000
        (50, 60]  0           0.833333
                  1           0.166667
        (60, 70]  0           0.666667
                  1           0.333333
3       (0, 5]    1           0.571429
                  0           0.428571
        (10, 15]  0           0.571429
                  1           0.428571
        (15, 20]  0           0.784615
                  1           0.215385
        (20, 25]  0           0.802817
                  1           0.197183
        (25, 30]  0           0.724138
                  1           0.275862
        (30, 40]  0           0.793651
                  1           0.206349
        (40, 50]  0           0.933333
                  1           0.066667
        (5, 10]   0           0.812500
                  1           0.187500
        (50, 60]  0           1.000000
        (60, 70]  0           0.666667
                  1           0.333333
        (70, 80]  0           1.000000
dtype: float64



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

Попробуем визуализацию от seaborn — уж очень красивые картинки она дает, хотя я больше привык к тексту.

sns.pairplot(train_data_munged, vars=["AgeF", "Pclass", "SexF"], hue="Survived", dropna=True)
sns.plt.show()




Красиво, но например корреляция в паре «класс-пол» не очень наглядна.

Оценим важность наших признаков алгоритмом SelectKBest.

selector = SelectKBest(f_classif, k=5)
selector.fit(train_data_munged[predictors], train_data_munged["Survived"])

scores = -np.log10(selector.pvalues_)

plt.bar(range(len(predictors)), scores)
plt.xticks(range(len(predictors)), predictors, rotation='vertical')
plt.show()




Вот вам статья с описанием того как именно он это делает. В параметрах SelectKBest можно указать и другие стратегии.

В принципе, мы и так уже все знаем — очень важен пол. Титулы важны — но у них сильная коррелляция с полом. Важен класс билета и каким-то образом палуба «F».

Оценка классификации



Перед тем как начинать запуск какой-либо классификации, нам нужно понять каким образом мы будем ее оценивать. В случае с конкурсами Kaggle все очень просто: мы просто читаем их правила. В случае Титаника оценкой будет служить отношение правильных оценок классификатора к общему числу пассажиров. Иными словами эта оценка называется accuracy.

Но прежде чем отправлять результат классификации по test-выборке на оценку в Kaggle, нам неплохо бы сначала самим понять хотя бы приблизительное качество работы нашего классификатора. Понять это мы сможем только используя train-выборку, поскольку только она содержит промаркированные данные. Но остается вопрос — как именно?

Часто в примерах можно увидеть нечто подобное:

classifier.fit(train_X, train_y)
predict_y = classifier.predict(train_X)
return metrics.accuracy_score(train_y, predict_y)


То есть мы обучаем классификатор на train-наборе, после чего на нем же его и проверяем. Несомненно в какой-то степени это дает некую оценку качества работы классификатора, но в целом этот подход неверен. Классификатор должен описывать не данные на котором его тренировали, а некую модель которая эти данные сгенерировала. В противном случае классификатор отлично подстраивается под train-выборку, при проверке на ней же показывает отличные результаты, однако при проверке на неком другом наборе данных с треском сливается. Что и называется overfitting.

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

Можно уже представить в голове циклы которые будут разделять данные, производить обучение и оценивание, но фишка в том, что все что нужно для реализации этого в sklearn — это определить стратегию.

cv = StratifiedKFold(train_data["Survived"], n_folds=3, shuffle=True, random_state=1)


Здесь мы определяем достаточно сложный процесс: тренировочные данные будут разделены на три куска, причем записи будут попадать в каждый кусок случайным образом (чтобы нивелировать возможную зависимость от порядка), кроме того стратегия отследит чтобы отношение классов в каждом куске было приблизительно равным. Таким образом мы будем производить три измерения по кускам 1+2 vs 3, 1+3 vs 2, 2+3 vs 1 — после этого сможем получить среднюю оценку аккуратности классификатора (что будет характеризовать качество работы), а также дисперсию оценки (что будет характеризовать стабильность его работы).

Классификация



Теперь протестируем работу различных классификаторов.

KNeighborsClassifier:
alg_ngbh = KNeighborsClassifier(n_neighbors=3)
scores = cross_val_score(alg_ngbh, train_data_scaled, train_data_munged["Survived"], cv=cv, n_jobs=-1)
print("Accuracy (k-neighbors): {}/{}".format(scores.mean(), scores.std()))


SGDClassifier:
alg_sgd = SGDClassifier(random_state=1)
scores = cross_val_score(alg_sgd, train_data_scaled, train_data_munged["Survived"], cv=cv, n_jobs=-1)
print("Accuracy (sgd): {}/{}".format(scores.mean(), scores.std()))


SVC:
alg_svm = SVC(C=1.0)
scores = cross_val_score(alg_svm, train_data_scaled, train_data_munged["Survived"], cv=cv, n_jobs=-1)
print("Accuracy (svm): {}/{}".format(scores.mean(), scores.std()))


GaussianNB:
alg_nbs = GaussianNB()
scores = cross_val_score(alg_nbs, train_data_scaled, train_data_munged["Survived"], cv=cv, n_jobs=-1)
print("Accuracy (naive bayes): {}/{}".format(scores.mean(), scores.std()))


LinearRegression:
def linear_scorer(estimator, x, y):
    scorer_predictions = estimator.predict(x)

    scorer_predictions[scorer_predictions > 0.5] = 1
    scorer_predictions[scorer_predictions <= 0.5] = 0

    return metrics.accuracy_score(y, scorer_predictions)

alg_lnr = LinearRegression()
scores = cross_val_score(alg_lnr, train_data_scaled, train_data_munged["Survived"], cv=cv, n_jobs=-1,
                         scoring=linear_scorer)
print("Accuracy (linear regression): {}/{}".format(scores.mean(), scores.std()))



Метод linear_scorer нужен поскольку LinearRegression — это регрессия возвращающая любое вещественное число. Соответственно мы разделяем шкалу границей 0.5 и приводим любые числа к двум классам — 0 и 1.

LogisticRegression:
alg_log = LogisticRegression(random_state=1)
scores = cross_val_score(alg_log, train_data_scaled, train_data_munged["Survived"], cv=cv, n_jobs=-1,
                         scoring=linear_scorer)
print("Accuracy (logistic regression): {}/{}".format(scores.mean(), scores.std()))


RandomForestClassifier:
alg_frst = RandomForestClassifier(random_state=1, n_estimators=500, min_samples_split=8, min_samples_leaf=2)
scores = cross_val_score(alg_frst, train_data_scaled, train_data_munged["Survived"], cv=cv, n_jobs=-1)
print("Accuracy (random forest): {}/{}".format(scores.mean(), scores.std()))


У меня получилось примерно так
Accuracy (k-neighbors): 0.698092031425/0.0111105442611
Accuracy (sgd): 0.708193041526/0.0178870678457
Accuracy (svm): 0.693602693603/0.018027360723
Accuracy (naive bayes): 0.791245791246/0.0244349506813
Accuracy (linear regression): 0.805836139169/0.00839878201296
Accuracy (logistic regression): 0.806958473625/0.0156323100754
Accuracy (random forest): 0.827160493827/0.0063488824349


Алгоритм Random Forest победил и дисперсия у него неплохая — кажется он стабилен.

Еще лучше



Вроде все хорошо и можно отправлять результат, но остался один мутный момент: у каждого классификатора есть свои параметры — как нам понять что мы выбрали наилучший вариант? Без сомнения можно долго сидеть и перебирать параметры вручную — но что если поручить эту работу компьютеру?

alg_frst_model = RandomForestClassifier(random_state=1)
alg_frst_params = [{
    "n_estimators": [350, 400, 450],
    "min_samples_split": [6, 8, 10],
    "min_samples_leaf": [1, 2, 4]
}]
alg_frst_grid = GridSearchCV(alg_frst_model, alg_frst_params, cv=cv, refit=True, verbose=1, n_jobs=-1)
alg_frst_grid.fit(train_data_scaled, train_data_munged["Survived"])
alg_frst_best = alg_frst_grid.best_estimator_
print("Accuracy (random forest auto): {} with params {}"
      .format(alg_frst_grid.best_score_, alg_frst_grid.best_params_))


Получается еще лучше!
Accuracy (random forest auto): 0.836139169473 with params {'min_samples_split': 6, 'n_estimators': 350, 'min_samples_leaf': 2}


Подбор можно сделать еще тоньше при наличии времени и желания — либо изменив параметры, либо используя другую стратегию подбора, например RandomizedSearchCV.

Пробуем xgboost



Все хвалят xgboost — давайте попробуем и его.

ald_xgb_model = xgb.XGBClassifier()
ald_xgb_params = [
    {"n_estimators": [230, 250, 270],
     "max_depth": [1, 2, 4],
     "learning_rate": [0.01, 0.02, 0.05]}
]
alg_xgb_grid = GridSearchCV(ald_xgb_model, ald_xgb_params, cv=cv, refit=True, verbose=1, n_jobs=1)
alg_xgb_grid.fit(train_data_scaled, train_data_munged["Survived"])
alg_xgb_best = alg_xgb_grid.best_estimator_
print("Accuracy (xgboost auto): {} with params {}"
      .format(alg_xgb_grid.best_score_, alg_xgb_grid.best_params_))


Почему-то тренировка зависала при использовании всех ядер, поэтому я ограничился одним потоком (n_jobs=1), но и в однопоточном режиме тренировка и классификация в xgboost работает очень быстро.

Результат тоже неплох
Accuracy (xgboost auto): 0.835016835017 with params {'n_estimators': 270, 'learning_rate': 0.02, 'max_depth': 2}


Результат



Классификатор выбран, параметры рассчитаны — осталось сформировать результат и отправить на проверку в Kaggle.

alg_test = alg_frst_best

alg_test.fit(train_data_scaled, train_data_munged["Survived"])

predictions = alg_test.predict(test_data_scaled)

submission = pd.DataFrame({
    "PassengerId": test_data["PassengerId"],
    "Survived": predictions
})

submission.to_csv("titanic-submission.csv", index=False)


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

  • Борьба в топе идет на сотые доли процентов — поэтому решает даже одна правильная или неправильная классификация из верификационного набора. Также на результат влияют случайные модификации параметров и алгоритма;
  • Лучший результат при локальной кросс-валидации не гарантирует лучшего результата при проверке верификационного набора. Бывает что реализация вроде бы разумной гипотезы улучшает локальный результат кросс-валидации и ухудшает результат проверки на Kaggle;
  • Вышеперечисленный два пункта приводят к третьему — скрипт и набор для отправки на верификацию должны находится под управлением любой системы управления версиями — в случае прорыва в счете нужно обязательно коммититься с указанием полученного результата в описании коммита;
  • Само соревнование несколько искусственно. В реальном мире у вас нет верификационного набора — только размеченная выборка на которой можно проводить кросс-валидацию, окончательное качество работы алгоритма оценивается как правило косвенно, например, по количеству жалоб пользователей;
  • Не очень понятна максимальная степень качества оценки которую можно получить машинным обучением на задаче — размер обучающего набора не очень высок, а сам процесс неочевиден. На самом деле — что такое процесс в случае бегства пассажиров с гибнущего корабля? Вполне вероятно, что не все усаженные в лодки выжили, также как вполне вероятно что не все неусаженные в лодки погибли — часть процесса должна быть очень хаотична и должна сопротивляться какому-либо моделированию.


Глава, в которой автор одновременно испытал удивление и просветление



Просматривая топ конкурсантов, нельзя не заметить людей которые набрали 1 (все ответы правильные) — причем у некоторых это получилось с первой же попытки.



На ум приходит следующий вариант: кто-то зарегистрировал учетку с которой начал перебором подбирать (позволяется не более 10 попыток в день) правильные ответы. Если я правильно понимаю это является разновидностью задач о взвешивании.

Однако подумав еще немного, нельзя не улыбнуться своей догадке: мы говорим о задаче, ответы по которой уже давно известны! В самом деле, гибель Титаника была потрясением для современников, и этому событию были посвящены фильмы, книги и документалистика. И скорее всего где-то есть полный поименный перечень пассажиров Титаника с описанием их судьбы. Но к машинному обучению это уже не относится.

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

Код



Полный код скрипта есть тут. Не забудьте выбрать признаки для обучения.

Комментарии (6)


  1. valodik
    29.12.2015 11:43

    Спасибо! Люблю, когда статья и сама по делу, и список близких хороших статей имеет.)


  1. Dumbris
    29.12.2015 12:02

    Заметил, что для random forest, в первом примере, количество n_estimators было задано 500. Grid search обнаружил более точный результат при n_estimators 350. Может попробовать варьировать n_estimators в более широком диапазоне?


    1. MzMz
      29.12.2015 13:29

      Первый вызов (500) — это просто иллюстрация того, как можно использовать классификатор, поэтому параметры взяты практически с потолка.

      В зачет Kaggle уходит классификатор из автоподбора GridSearchCV — я там менял диапазоны и в более узких и в более широких параметрах, правда без особого фанатизма, это быстро надоедает :)


  1. ternaus
    29.12.2015 16:01

    Замечательно пишете.

    Вы в других соревнованиях не участвовали? Титаник это хорошо, но просто. Было бы интересно почитать ваш разбор какого-нибудь более насыщенного соревнования.


  1. ciiccii
    29.12.2015 19:23

    Пример для xgboost у меня почему-то не работает. Падает с ошибкой:

    XGBoostError: b'base_score must be in (0,1) for logistic loss'
    
    на строке alg_xgb_grid.fit(train_data_scaled, train_data_munged[«Survived»])

    В чём может быть ошибка?


    1. MzMz
      29.12.2015 20:37

      github.com/dmlc/xgboost/blob/master/src/learner/objective-inl.hpp#L104

      могу только предположить, что либо train_data_munged[«Survived»] вышел за пределы, либо у классификатора есть какой-то незаданный параметр который в вашей версии библиотеки (у меня 0.4) установлен в другое значение.