{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "## Открытый курс по машинному обучению. Сессия № 3\n", "\n", "###
Автор материала: Михаил Каменев " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##
Индивидуальный проект по анализу данных
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**План исследования**\n", " - Описание набора данных и признаков\n", " - Первичный анализ признаков\n", " - Первичный визуальный анализ признаков\n", " - Закономерности, \"инсайты\", особенности данных\n", " - Предобработка данных\n", " - Создание новых признаков и описание этого процесса\n", " - Кросс-валидация, подбор параметров\n", " - Построение кривых валидации и обучения \n", " - Прогноз для тестовой или отложенной выборки\n", " - Оценка модели с описанием выбранной метрики\n", " - Выводы\n", " \n", " Более детальное описание [тут](https://goo.gl/cJbw7V)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "from matplotlib import pyplot as plt\n", "import seaborn as sns\n", "import scipy\n", "from statsmodels.stats.weightstats import *\n", "\n", "from sklearn.linear_model import RidgeCV, Ridge, Lasso, LassoCV\n", "from sklearn.ensemble import RandomForestRegressor\n", "from sklearn.model_selection import train_test_split, KFold\n", "from sklearn.preprocessing import StandardScaler, LabelBinarizer, PolynomialFeatures\n", "from sklearn.metrics import mean_absolute_error, make_scorer\n", "from xgboost import XGBClassifier\n", "from hyperopt import fmin,tpe, hp, STATUS_OK, Trials\n", "import xgboost as xgb\n", "from sklearn.model_selection import learning_curve, validation_curve\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Часть 1. Описание набора данных и признаков" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "Датасет содержит информацию о 53940 бриллиантах. По некоторым характеристикам (об этом позже) будем предсказывать стоимость. Данные можно скачать здесь.\n", "\n", "С точки зрения бизнеса, ценность задачи понятна - по характеристикам бриллианта предсказать, сколько долларов за него можно получить. От бизнеса я далёк, поэтому интерес чисто спортивный: разобраться, какие характеристики и как влияют на стоимость этих камешков =)\n", "\n", "Признаки\n", "- carat - вес бриллианта в каратах, вещественный\n", "- cut - качество огранки, категориальный. Принимает пять возможных значений: Fair, Good, Very Good, Premium, Ideal\n", "- color - \"цвет\" бриллианта. Категориальный признак, принимает значения J,I,H,G,F,E,D (от худшего (J) к лучшему (D))\n", "- clarity - чистота бриллианта. Категориальный признак, принимает значения I1 (худший), SI2, SI1, VS2, VS1, VVS2, VVS1, IF (лучший)\n", "- x,y,z - три признака, характеризущие размеры бриллианта, вещественные\n", "- depth - признак, который высчитывается на основе трех предыдущих по формуле 2 * z / (x + y), вещественный\n", "- table - отношение ширины верхней грани бриллианты к его максимальной ширине, в процентах\n", "\n", "\n", "Целевой признак: price - стоимость бриллианта в долларах\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Часть 2. Первичный анализ признаков" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#загружаем dataset\n", "diamonds_df = pd.read_csv('../../data/diamonds.csv')\n", "diamonds_df.head()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "diamonds_df.describe()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Видно, что масштабы признаков отличаются. В дальнейшем нужно будет применить StandartScaler" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "diamonds_df.info()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "В данных отсутствуют пропуски. Итого, имеется 6 вещественных, 1 целочисленный (unnamed: 0 не считаем) и 3 категориальных признака." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Анализ целочисленных и вещественных признаков" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "real_features = ['carat', 'depth', 'table', 'x', 'y', 'z','price']" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Изучим корреляцию вещественных признаков и целевой переменной\n", "sns.heatmap(diamonds_df[real_features].corr(method='spearman'));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Признаки carat, x,y,z имеют большую корреляцию, как между собой, так и с целевой переменной, что не удивительно. При этом, корреляция целевой переменной и признаков depth, table почти отсутствует" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Анализ категориальных признаков" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cat_features = ['cut','color','clarity']\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(16, 10))\n", "\n", "for idx, feature in enumerate(cat_features):\n", " sns.countplot(diamonds_df[feature], ax=axes[idx % 3], label=feature)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Реальные значения категориальных признаков не отличаются от тех, что заявлены в описании. Кроме того, видно, что уникальных значений не много, так что One Hot encoding должен отлично сработать. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Анализ целевого признака" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sns.distplot(diamonds_df['price'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Распределение имеет тяжелый правый хвост. Применим логарифмирование." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sns.distplot(diamonds_df['price'].map(np.log1p))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Помогло это не сильно: получилось бимодальное распределение. Зато хвост исчез =) Для наглядности, построим QQ график" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "stats.probplot(diamonds_df['price'], dist=\"norm\", plot=plt);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Выводы\n", "- Вещественные признаки (carat, depth, table, x, y, z) масштабируем\n", "- К категориальным признакам ('cut','color','clarity') применяем one hot encoding\n", "- Целевую переменную логарифмируем" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Часть 3. Первичный визуальный анализ признаков" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "#### Анализ целочисленных и вещественных признаков" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Начнем с построения гистограмм вещественных признаков\n", "fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(16, 10))\n", "\n", "for idx, feature in enumerate(real_features[:-1]): #price рисовать не будем\n", " sns.distplot(diamonds_df[feature], ax=axes[idx // 3, idx % 3], label=feature)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Распределение признаков depth, table, y, z отдаленно, но напоминает колокол. У depth хвосты тяжеловаты для нормального распределения; carat и table скорее бимодальные. Кроме того, у них тяжелые правые хвосты, так что np.log1p не помешает. По графикам выше не видно выбросов. Проверим, что это действительно так, с помощью boxplot" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(16, 10))\n", "\n", "for idx, feature in enumerate(real_features[:-1]): #price рисовать не будем\n", " sns.boxplot(diamonds_df[feature], ax=axes[idx // 3, idx % 3], orient='v')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Каких-либо серьезных аномалий в рассматриваемых данных нет. На всякий случай посмотрим бриллиант с y=60, z = 32 и carat > 4. Если он стоит дорого, то за выброс его считать не будем." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "diamonds_df[diamonds_df['y'] > 55].head()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "diamonds_df[diamonds_df['z'] > 30].head()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "diamonds_df[diamonds_df['carat'] > 4].head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Видно, что это просто очень дорогие камни. Посмотрим, как рассматриваемые признаки взаимосвязаны с целевой перменной" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sns.pairplot(diamonds_df[real_features], diag_kind=\"kde\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- вес бриллианта показывает степенную зависимость от его размеров\n", "- depth и table почти никак не взаимосвязаны с остальными признаками, в том числе и целевым\n", "- x,y,z связаны между собой линейно\n", "- цена линейно зависит от размеров\n", "- зависимость между ценой и весом сложно назвать линейной, но монотонный тренд есть" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Анализ категориальных признаков" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Посмотрим, как целевая переменная зависит от категориальных признаков" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# цвет бриллианта\n", "fig, axes = plt.subplots(nrows=3, ncols=3, figsize=(16, 10))\n", "\n", "for idx, (color, sub_df) in enumerate(pd.groupby(diamonds_df, 'color')): \n", " ax = sns.distplot(sub_df['price'], kde=False, ax=axes[idx // 3, idx % 3])\n", " ax.set(xlabel=color)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Распределения для всех значений цветов имеют тяжелый правый хвост и не сильно отличаются друг от друга." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# чистота бриллианта\n", "fig, axes = plt.subplots(nrows=3, ncols=3, figsize=(16, 10))\n", "\n", "for idx, (clarity, sub_df) in enumerate(pd.groupby(diamonds_df, 'clarity')): \n", " ax = sns.distplot(sub_df['price'], kde=False, ax=axes[idx // 3, idx % 3])\n", " ax.set(xlabel=clarity)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Хвосты у всех тяжелые, но у SI1,SI2 присутствуют дополнительные пики в районе 5000." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# качество огранки\n", "fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(16, 10))\n", "\n", "for idx, (cut, sub_df) in enumerate(pd.groupby(diamonds_df, 'cut')): \n", " ax = sns.distplot(sub_df['price'], kde=False, ax=axes[idx // 3, idx % 3])\n", " ax.set(xlabel=cut)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "И снова пики в районе 5000 (у Good и Premium). А в целом графики похожи." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Нарисуем boxplot для каждого значения" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(nrows=3, ncols=1, figsize=(16, 10))\n", "\n", "# Отобразим строки в числа в порядке от худшего к лучшему. Так удобнее на графике смотреть\n", "df = diamonds_df.copy()\n", "df['color'] = df['color'].map({'J': 0, 'I': 1, 'H': 2, 'G': 3, 'F': 4, 'E': 5, 'D': 6})\n", "df['clarity'] = df['clarity'].map({'I1': 0, 'SI2': 1, 'SI1': 2, 'VS2': 3, 'VS1': 4, 'VVS2': 5, 'VVS1': 6, 'IF': 7 })\n", "df['cut'] = df['cut'].map({'Fair': 0, 'Good': 1, 'Very Good': 2, 'Premium': 3, 'Ideal': 4})\n", "\n", "for idx, feature in enumerate(cat_features):\n", " sns.boxplot(x=feature, y='price',data=df,hue=feature, ax=axes[idx])\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Тут уже интереснее. Начнем с огранки. Видно, что медиана максимальна для Very Good и Premium. Для ideal медианное значение цены гораздо меньше. Аналогичные наблюдения можно сделать для цвета и чистоты. Возможно, бриллианты с наилучшими свойствами на очень большие, и, соответсвенно, их цена ниже. Проверим это." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(nrows=3, ncols=1, figsize=(16, 10))\n", "\n", "for idx, feature in enumerate(cat_features):\n", " sns.boxplot(x=feature, y='carat',data=df,hue=feature, ax=axes[idx])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Действительно, медианное значение веса для бриллиантов с очень хорошими характеристиками меньше, чем для бриллиантов с плохими харакетристиками. Напоследок, посмотрим сколько бриллиантов с той или иной харакеристикой присутствует в данных." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(16, 10))\n", "\n", "for idx, feature in enumerate(cat_features):\n", " sns.countplot(df[feature], ax=axes[idx % 3], label=feature)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Видно, что очень мало камней с плохой огранкой. Также мало камней с плохими цветовыми харакеристиками. Но и не очень много с идеальными. Распределение чистоты камня напоминает лапласовское распределение." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Часть 4. Закономерности, \"инсайты\", особенности данных" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "#### Основные выводы по предыдущим пунктам:\n", "- Ключевые признаки для прогнозирования: вес и размеры бриллианта (carat, x, y, z). По графикам видно, что есть монотонная зависимость этих признаков и цены. Что логично\n", "- Признаки depth и table почти не влияют на стоимость камня\n", "- Исключительно по категориальным признакам сложно что-либо сказать о целевой переменной. Однако видно, что чем лучше бриллиант с точки зрения этих признаков, тем больше вероятность того, что он будет не очень большого размера\n", "- Выбросы в данных отсутствуют\n", "- Так как у целевой переменной очень тяжелый правый хвост, в качестве метрики будем использовать среднюю абсолютную ошибку, а не квадратичную.\n", "- Видно, что зависимость от ключевых признаков близка к линейной. Поэтому в качестве бейзлайна будем использовать линейную регрессию.\n", "- Более того, признаков не так уж и много, поэтому будем рассматривать также случайный лес и градиентный бустинг (тут он должен затащить =)). А случайный лес инетересен исключительно для сравнения с бустингом" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Часть 5. Предобработка данных " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Для начала, выделим выборку для тестирования\n", "X = diamonds_df.drop(['price'], axis=1).values[:,1:] # отсекаем индекс\n", "y = diamonds_df['price']\n", "\n", "X_train, X_test, y_train, y_test = train_test_split(X,y,test_size=0.3,random_state=4444, shuffle=True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# признаки с индексами 1,2,3 категориальные. Применим к ним ohe\n", "label_bin = LabelBinarizer()\n", "X_train_cut_ohe = label_bin.fit_transform(X_train[:,1])\n", "X_test_cut_ohe = label_bin.transform(X_test[:,1])\n", " \n", "X_train_color_ohe = label_bin.fit_transform(X_train[:,2])\n", "X_test_color_ohe = label_bin.transform(X_test[:,2])\n", " \n", "X_train_clarity_ohe = label_bin.fit_transform(X_train[:,3])\n", "X_test_clarity_ohe = label_bin.transform(X_test[:,3]) \n", "\n", "# carat, x и целевую переменную логарифмируем\n", "log_vect = np.vectorize(np.log1p)\n", "X_train_сarat_log = log_vect(X_train[:,0]).reshape(-1,1)\n", "X_test_сarat_log = log_vect(X_test[:,0]).reshape(-1,1)\n", "X_train_x_log = log_vect(X_train[:,6]).reshape(-1,1)\n", "X_test_x_log = log_vect(X_test[:,6]).reshape(-1,1)\n", "y_train_log = log_vect(y_train)\n", "y_test_log = log_vect(y_test)\n", "\n", "# масштабириуем вещественные признаки\n", "scaler = StandardScaler()\n", "X_train_real = np.hstack((X_train_сarat_log, X_train_x_log, X_train[:,[7,8,4,5]]))\n", "X_test_real = np.hstack((X_test_сarat_log, X_test_x_log, X_test[:,[7,8,4,5]]))\n", "X_train_real_scaled = scaler.fit_transform(X_train_real)\n", "X_test_real_scaled = scaler.transform(X_test_real)\n", "\n", "# В качестве дополнительных признаков будем рассматривать полиномиальные признаки\n", "#Данные признаки должны улучшить качество линейной модели.\n", "\n", "X_train_additional = PolynomialFeatures().fit_transform(X_train_real)\n", "X_test_additional = PolynomialFeatures().fit_transform(X_test_real)\n", "X_train_additional_scaled = scaler.fit_transform(X_train_additional)\n", "X_test_additional_scaled = scaler.transform(X_test_additional)\n", "\n", "\n", "\n", "# Объединяем все преобразованные признаки\n", "X_train_transformed = np.hstack((X_train_real_scaled,X_train_cut_ohe, X_train_color_ohe, X_train_clarity_ohe))\n", "X_test_transformed = np.hstack((X_test_real_scaled,X_test_cut_ohe, X_test_color_ohe, X_test_clarity_ohe))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Часть 6. Создание новых признаков и описание этого процесса" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "Смотри предыдущий пункт" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Часть 7. Кросс-валидация, подбор параметров" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "Рассмотрим сначала линейную модель. Данные разделим на 5 фолдов. С помощью RidgeCV и LassoCV будем оптимизировать силу регуляризации." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# функция потерь для рассматриваемой задачи. Ошибку смотрим на исходных данных\n", "def mean_absolute_exp_error(model, X,y):\n", " return -mean_absolute_error(np.expm1(model.predict(X)), np.expm1(y))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cv = KFold(n_splits=5, shuffle=True, random_state=4444)\n", "alphas = np.logspace(-5,2,100)\n", "ridge_cv = RidgeCV(alphas=alphas, scoring=mean_absolute_exp_error, cv=cv)\n", "lasso_cv = LassoCV(alphas=alphas, cv=cv, random_state=4444)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ridge_cv.fit(X_train_transformed, y_train_log)\n", "lasso_cv.fit(X_train_transformed, y_train_log)\n", "print('Optimized alpha: Ridge = %f, Lasso = %f' % (ridge_cv.alpha_, lasso_cv.alpha_))\n", "score_ridge = mean_absolute_error(y_test, np.expm1(ridge_cv.predict(X_test_transformed)))\n", "score_lasso = mean_absolute_error(y_test, np.expm1(lasso_cv.predict(X_test_transformed)))\n", "print('Ridge regression score = %f' % score_ridge)\n", "print('Lasso regression score = %f' % score_lasso)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Оба метода показали схожий результат. Что будет, если мы добавим новые признаки?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "X_train_transformed_add = np.hstack((X_train_transformed, X_train_additional_scaled))\n", "X_test_transformed_add = np.hstack((X_test_transformed, X_test_additional_scaled))\n", "ridge_cv.fit(X_train_transformed_add, y_train_log)\n", "lasso_cv.fit(X_train_transformed_add, y_train_log)\n", "print('Optimized alpha: Ridge = %f, Lasso = %f' % (ridge_cv.alpha_, lasso_cv.alpha_))\n", "score_ridge = mean_absolute_error(y_test, np.expm1(ridge_cv.predict(X_test_transformed_add)))\n", "score_lasso = mean_absolute_error(y_test, np.expm1(lasso_cv.predict(X_test_transformed_add)))\n", "print('Ridge regression score = %f' % score_ridge)\n", "print('Lasso regression score = %f' % score_lasso)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ошибка значительно уменьшилась. Построим кривые валидации и обучения" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "# код из статьи на хабре\n", "\n", "\n", "def plot_with_err(x, data, **kwargs):\n", " mu, std = data.mean(1), data.std(1)\n", " lines = plt.plot(x, mu, '-', **kwargs)\n", " plt.fill_between(x, mu - std, mu + std, edgecolor='none',\n", " facecolor=lines[0].get_color(), alpha=0.2)\n", "\n", "model = Ridge(random_state=4444) \n", "alphas = np.logspace(1,2,10) + 10 # Если коэффициент регуляризации мал, то значения получаются заоблочными\n", "val_train, val_test = validation_curve(model, X_train_transformed_add, y_train_log,'alpha', alphas, cv=cv,scoring=mean_absolute_exp_error)\n", "plot_with_err(alphas, -val_train, label='training scores')\n", "plot_with_err(alphas, -val_test, label='validation scores')\n", "plt.xlabel(r'$\\alpha$'); plt.ylabel('MAE')\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Судя по кривым валидации, модель недообучилась: ошибки лежат близко друг к другу. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# код из статьи на хабре\n", " \n", "def plot_learning_curve(model, X,y):\n", " train_sizes = np.linspace(0.05, 1, 20)\n", " \n", " N_train, val_train, val_test = learning_curve(model,X, y, train_sizes=train_sizes, cv=5,scoring=mean_absolute_exp_error, random_state=4444)\n", " plot_with_err(N_train, -val_train, label='training scores')\n", " plot_with_err(N_train, -val_test, label='validation scores')\n", " plt.xlabel('Training Set Size'); plt.ylabel('MAE')\n", " plt.legend()\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model = Ridge(alpha=52.140083,random_state=4444)\n", "plot_learning_curve(model, X_train_transformed_add, y_train_log)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Кривые лежат близко друг к другу почти с самого начала. Вывод: наблюдений у нас достаточно, нужно двигаться в сторону усложнения модели" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Случайный лес" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Случайный лес должен хорошо работать \"из коробки\". Поэтому будем оптимизировать только число деревьев.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "model = RandomForestRegressor(n_estimators=100, random_state=4444)\n", "n_estimators = [10,25,50,100,250,500,1000]\n", "val_train, val_test = validation_curve(model, X_train_transformed, y_train_log,'n_estimators', n_estimators, cv=cv,scoring=mean_absolute_exp_error)\n", " \n", "plot_with_err(n_estimators, -val_train, label='training scores')\n", "plot_with_err(n_estimators, -val_test, label='validation scores')\n", "plt.xlabel('n_estimators'); plt.ylabel('MAE')\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Видно, что начиная с 200 деревьев качество практически не изменяется. Поэтому в качестве еще одной модели будем рассматривать случайный лес именно с таким количеством деревьев." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "forest_model = RandomForestRegressor(n_estimators=200, random_state=4444)\n", "forest_model.fit(X_train_transformed, y_train_log)\n", "forest_prediction = np.expm1(forest_model.predict(X_test_transformed))\n", "score = mean_absolute_error(y_test, forest_prediction)\n", "print('Random forest score: %f' % score)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# посмотрим на важность признаков\n", "np.argsort(forest_model.feature_importances_)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Первые четыре столбца обучающей выборки соответствуют признакам carat, x,y,z. Как и предполагалось в начале, 3 из 4 этих признаков имеют наибольшую важность для модели" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "# Построим, также, кривую обучения\n", "plot_learning_curve(model, X_train_transformed, y_train_log)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "График выходит на полку, так что больше данных нам не нужно" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### boosting. А что boosting?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "X_train_boosting, X_valid_boosting, y_train_boosting, y_valid_boosting = train_test_split(\n", " X_train_transformed, y_train_log, test_size=0.3, random_state=4444)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def score(params):\n", " from sklearn.metrics import log_loss\n", " print(\"Training with params:\")\n", " print(params)\n", " params['max_depth'] = int(params['max_depth'])\n", " dtrain = xgb.DMatrix(X_train_boosting, label=y_train_boosting)\n", " dvalid = xgb.DMatrix(X_valid_boosting, label=y_valid_boosting)\n", " model = xgb.train(params, dtrain, params['num_round'])\n", " predictions = model.predict(dvalid).reshape((X_valid_boosting.shape[0], 1))\n", " score = mean_absolute_error(np.expm1(y_valid_boosting), np.expm1(predictions))\n", " # score = mean_absolute_error(y_valid_boosting, predictions)\n", " print(\"\\tScore {0}\\n\\n\".format(score))\n", " return {'loss': score, 'status': STATUS_OK}\n", "\n", "def optimize(trials):\n", " space = {\n", " 'num_round': 200,\n", " 'learning_rate': hp.quniform('eta', 0.05, 0.5, 0.005),\n", " 'max_depth': hp.quniform('max_depth', 3, 14, 1),\n", " 'min_child_weight': hp.quniform('min_child_weight', 1, 10, 1),\n", " 'subsample': hp.quniform('subsample', 0.5, 1, 0.05),\n", " 'gamma': hp.quniform('gamma', 0.5, 1, 0.01),\n", " 'colsample_bytree': hp.quniform('colsample_bytree', 0.4, 1, 0.05),\n", " 'eval_metric': 'mae',\n", " 'objective': 'reg:linear',\n", " 'nthread' : 4,\n", " 'silent' : 1,\n", " 'seed': 4444\n", " }\n", " \n", " best = fmin(score, space, algo=tpe.suggest,trials=trials, max_evals=100)\n", " return best" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "# Оптимизация параметров\n", "trials = Trials()\n", "best_params = optimize(trials)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "best_params" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "params = {\n", " 'num_round': 200,\n", " 'colsample_bytree': 0.65,\n", " 'eta': 0.145,\n", " 'gamma': 0.55,\n", " 'max_depth': 10,\n", " 'min_child_weight': 4.0,\n", " 'subsample': 1.0,\n", " 'eval_metric': 'mae',\n", " 'objective': 'reg:linear',\n", " 'nthread' : 4,\n", " 'silent' : 1,\n", " 'seed': 4444}" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dtrain = xgb.DMatrix(X_train_transformed, label=y_train_log)\n", "dvalid = xgb.DMatrix(X_test_transformed, label=y_test_log)\n", "boosting_model = xgb.train(params, dtrain, params['num_round'])\n", "predictions = boosting_model.predict(dvalid).reshape((X_test_transformed.shape[0], 1))\n", "score = mean_absolute_error(y_test, np.expm1(predictions))\n", "print('Boosting score: %f' % score)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Часть 8. Построение кривых валидации и обучения " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "В большом количестве в предыдущем пункте" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Часть 9. Прогноз для тестовой или отложенной выборки" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "В большом количестве в части 7" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Часть 10. Оценка модели с описанием выбранной метрики" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "Приведем результаты различных моделей на тестовой выборке.\n", "Как уже оговаривалось ранее, в качестве метрики используем MAE" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "pure_ridge = Ridge(random_state=4444, alpha=0.00001) # гребневая регрессия на исходных данных\n", "pure_ridge.fit(X_train_transformed, y_train_log)\n", "pure_ridge_score = mean_absolute_error(y_test, np.expm1(pure_ridge.predict(X_test_transformed)))\n", "print('Ridge regression score: %f' % pure_ridge_score)\n", "poly_ridge = Ridge(random_state=4444, alpha=52.140083) # гребневая регрессия с полиномиальными признаками\n", "poly_ridge.fit(X_train_transformed_add, y_train_log)\n", "poly_ridge_score = mean_absolute_error(y_test, np.expm1(poly_ridge.predict(X_test_transformed_add)))\n", "print('Ridge regression score with poly features: %f' % poly_ridge_score)\n", "forest_score = mean_absolute_error(y_test, np.expm1(forest_model.predict(X_test_transformed)))\n", "print('Random forest score: %f' % forest_score)\n", "boosting_score = mean_absolute_error(y_test, np.expm1(boosting_model.predict(dvalid)))\n", "print('XGBoost score: %f' % boosting_score)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Результаты близки к тем, что получались на кросс-валидации. Так что всё хорошо =)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Часть 11. Выводы " ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "В данном проекте рассматривались достаточно \"простые\" данные, поэтому основной упор был сделан на применение различных моделей для их анализа. С одной стороны, случайный лес без какой-либо настройки гиперпараметров показал лучший результат. С другой стороны, если потратить больше времени на оптимизацию градиентного бустинга, возможно, он сможет показать результат лучше, чем у случайного леса. Стоит, также, отметить линейную модель: после добавления полиномиальных признаков она показала очень неплохой результат (если сравнивать с моделью без дополнительных признаков =)). Зато сложность гораздо меньше. Если вдруг кому-то по жизни придется оценивать бриллианты, можете смело использовать предложенную модель случайного леса. В среднем будете терять по 275 $ с одного камушка :p\n", "\n", "Спасибо за внимание!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.4" } }, "nbformat": 4, "nbformat_minor": 1 }