В общем эта статья для еще начинающих от уже начавшего.
Титаник
Включаем подходящую для работы музыку и начинаем исследование.
Конкурс «по титанику» был уже неоднократно отмечен на Хабре. Особенно хотелось бы отметить последнюю статью из списка — оказывается, что исследование данных может быть не менее интригующим чем хороший детективный роман, а отличный результат (0.81340) можно получить не только Random Forest классификатором.
- habrahabr.ru/company/microsoft/blog/268039 — Предсказание выживания пассажиров Титаника при помощи Azure Machine Learning
- habrahabr.ru/post/165001 — Data Mining: Первичная обработка данных при помощи СУБД
- habrahabr.ru/company/mlclass/blog/248779 — Когда данных действительно много: Vowpal Wabbit
- habrahabr.ru/post/272201 — Устойчивая красота неприличных моделей
- habrahabr.ru/post/202090 — Основы анализа данных на python с использованием pandas+sklearn
- habrahabr.ru/company/mlclass/blog/270973 — Титаник на Kaggle: вы не дочитаете этот пост до конца
Также хотелось бы отметить статью о другом конкурсе. Из нее можно понять каким именно образом должен работать мозг исследователя и что большая часть времени должна быть уделена предварительному анализу и обработке данных.
- habrahabr.ru/post/270367 — Как я победил в конкурсе BigData от Beeline
Инструментарий
Для решения задачи я использую 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ам я их не пробовал, но на первый взгляд выглядят они многообещающе.
- www.continuum.io/downloads — Anaconda
- www.enthought.com/products/canopy — Canopy
Данные
Скачаем исходные данные для задачи и посмотрим что же нам выдали.
$ 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
Видим следующее:
- Возраст заполнен не у всех
- Билеты имеют какой-то странный и неконсистентный формат
- В именах есть титул (мисс, мистер, миссис и т.д.)
- Номера кают прописаны мало у кого (тут есть леденящая душу история почему)
- В существующих номерах кают видимо прописан код палубы (так и есть, как оказалось)
- Также, согласно статье, в номере каюты зашифрована сторона
- Сортируем по имени. Видно, что многие путешествовали семьями, и виден масштаб трагедии — часто семьи оказались разделены, выжила лишь часть
- Сортируем по билету. Видно, что по одному коду билета путешествовало сразу несколько людей, причем часто — с разными фамилиями. Беглый взгляд вроде бы показывает, что люди с одинаковым номером билета часто разделяют одну и ту же участь.
- У части пассажиров не проставлен порт посадки
Вроде приблизительно все ясно, переходим непостредственно к работе с данными.
Загрузка данных
Чтобы не зашумлять дальшейший код сразу приведу все используемые импорты:
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 (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_))
Подбор можно сделать еще тоньше при наличии времени и желания — либо изменив параметры, либо используя другую стратегию подбора, например 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 работает очень быстро.
Результат
Классификатор выбран, параметры рассчитаны — осталось сформировать результат и отправить на проверку в 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)
Dumbris
29.12.2015 12:02Заметил, что для random forest, в первом примере, количество n_estimators было задано 500. Grid search обнаружил более точный результат при n_estimators 350. Может попробовать варьировать n_estimators в более широком диапазоне?
MzMz
29.12.2015 13:29Первый вызов (500) — это просто иллюстрация того, как можно использовать классификатор, поэтому параметры взяты практически с потолка.
В зачет Kaggle уходит классификатор из автоподбора GridSearchCV — я там менял диапазоны и в более узких и в более широких параметрах, правда без особого фанатизма, это быстро надоедает :)
ternaus
29.12.2015 16:01Замечательно пишете.
Вы в других соревнованиях не участвовали? Титаник это хорошо, но просто. Было бы интересно почитать ваш разбор какого-нибудь более насыщенного соревнования.
ciiccii
29.12.2015 19:23Пример для xgboost у меня почему-то не работает. Падает с ошибкой:
на строке alg_xgb_grid.fit(train_data_scaled, train_data_munged[«Survived»])XGBoostError: b'base_score must be in (0,1) for logistic loss'
В чём может быть ошибка?MzMz
29.12.2015 20:37github.com/dmlc/xgboost/blob/master/src/learner/objective-inl.hpp#L104
могу только предположить, что либо train_data_munged[«Survived»] вышел за пределы, либо у классификатора есть какой-то незаданный параметр который в вашей версии библиотеки (у меня 0.4) установлен в другое значение.
valodik
Спасибо! Люблю, когда статья и сама по делу, и список близких хороших статей имеет.)