{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Открытый курс по машинному обучению. Сессия № 3\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Автор материала: Романова Анастасия (Slack: @anastasiaromane)\n",
    "#### Прогнозирование посещаемости в национальных парках и др. достопримечательностях США "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**План исследования**\n",
    " - Часть 1. Описание набора данных и признаков\n",
    " - Часть 2-3. Первичный и визуальный анализ признаков\n",
    " - Часть 4. Закономерности, \"инсайты\", особенности данных\n",
    " - Часть 5-6. Выбор метрики и модели\n",
    " - Часть 7-8. Предобработка данных и создание новых признаков и описание этого процесса\n",
    " - Часть 9. Кросс-валидация, подбор параметров\n",
    " - Часть 10. Построение кривых валидации и обучения \n",
    " - Часть 11. Прогноз для тестовой или отложенной выборки. Оценка модели с описанием выбранной метрики\n",
    " - Часть 12. Выводы"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Импорт библиотек"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.stats as stats\n",
    "from matplotlib import pyplot as plt\n",
    "import seaborn as sns\n",
    "from sklearn.preprocessing import OneHotEncoder, LabelEncoder\n",
    "from sklearn.linear_model import Ridge\n",
    "from sklearn.model_selection import train_test_split, KFold, GridSearchCV, validation_curve, learning_curve\n",
    "from catboost import Pool, CatBoostRegressor\n",
    "from sklearn.metrics import r2_score\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 1. Описание набора данных и признаков"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Задача"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Задача данного проекта - предсказать количество посетителей в национальных парках и других достопримечательностях США для последнего доступного года в датасете. Будем рещать задачу регрессии.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Данные"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Данные взяты с этого [сайта](https://data.world/inform8n/us-national-parks-visitation-1904-2016-with-boundaries). Датасет включает в себя исторические данные о посещаемости национальных парков и достопримечательной США с 1904 по 2016 год."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Признаки"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Загрузим наш основной датасет."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PATH_TO_DATA = ('data')\n",
    "df = pd.read_csv(os.path.join(PATH_TO_DATA, 'national-parks-visitation.csv'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.head(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Полный датасет состоит из 21560 объектов и 18 признаков."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- **Created By** - метаданные датасета\n",
    "- **Measure Selector** - метаданные датасета\n",
    "- **Year** - год в формате даты\n",
    "- **Date Edit** - метаданные датасета\n",
    "- **ScrapeURL** - метаданные датасета\n",
    "- **GIS Notes** - метаданные датасета\n",
    "- **Gnis Id** - id достопримечатльности из базы данных Geographic Names Information System \n",
    "- **Geometry** - геометрическая форма парка\n",
    "- **Metadata** - метаданные датасета\n",
    "- **Number of Records** - метаданные датасета\n",
    "- **Parkname** - название достопримечатльности (colloquial)\n",
    "- **Region** - район в системе национальных парков США (подробнее [здесь](https://en.wikipedia.org/wiki/Organization_of_the_National_Park_Service )) \n",
    "- **State** - штат, где расположена достопримечатльность\n",
    "- **Unit Code** - уникальный код достопримечатльности\n",
    "- **Unit Name** - полное название достопримечатльности\n",
    "- **Unit Type** - тип достопримечатльности\n",
    "- **Visitors** - метаданные датасета (**Целевой признак**)\n",
    "- **YearRaw** - год в формате текста"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Загрузим дополнительный датасет, чтобы по GNIS ID вытащить широту и долготу объектов. Данный датасет взят с этого [сайта](https://geonames.usgs.gov/domestic/download_data.htm)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gnis_df = pd.read_csv(os.path.join(PATH_TO_DATA, 'gnis-april2018.txt'), sep='|').rename(columns={'FEATURE_ID':'Gnis Id'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gnis_df.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Удалим ненужные колонки, преобразуем GNIS ID в текстовый формат и объединим два датасета по GNIS ID."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gnis_df.drop(['FEATURE_NAME', 'FEATURE_CLASS', 'STATE_ALPHA','STATE_NUMERIC',\n",
    "              'COUNTY_NAME', 'COUNTY_NUMERIC', 'PRIMARY_LAT_DMS','PRIM_LONG_DMS',\n",
    "              'SOURCE_LAT_DMS', 'SOURCE_LONG_DMS', 'SOURCE_LAT_DEC','SOURCE_LONG_DEC',\n",
    "              'ELEV_IN_M','ELEV_IN_FT', 'MAP_NAME', 'DATE_CREATED','DATE_EDITED'], axis=1, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gnis_df['Gnis Id'] = gnis_df['Gnis Id'].astype('str')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.merge(df, gnis_df, how='left', on='Gnis Id')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 2-3. Первичный + визуальный анализ признаков"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим на тип признаков в главном датасете."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Из 18 признаков только 4 вещественных, включая целевой признак Visitors, все остальные категориальные или текстовые.\n",
    "Посмотрим еще раз на данные."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Первым делом необходимо удалить некоторые признаки из датасета, так как большинство из них являются метаданными датасета и не несут большой информативности для анализа. Кроме того, удалим признак Year в формате даты, так как нам будет достатотчно признака YearRaw."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.drop(['Created By', 'Measure Selector', 'Year','Date Edit', 'ScrapeURL',\n",
    "       'GIS Notes', 'Geometry', 'Metadata', 'Number of Records'], axis=1, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Проверим данные на пропуски."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.isnull().sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4 признака имеют пропуски. Для каждого из признаков разные решения. Например, для признака Parkname можно заменить пропуски на Unknown. На самом деле, этот признак можно удалить (Сделаем это чуть позднее), т.к. есть другой более точной признак, который содержит в себе полное название достопримечательни Unit Name и не содержит пропусков."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['Parkname'].fillna('Unknown', inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим на другие пропуски в датасете."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[df.isnull().any(axis=1)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "У двух объектов отсутсвует широта и долгота. Проблема в том, что эти два объекта не имеют GNIS ID, по которому мы находили географические данные. В таком случае можно поискать долготу и широту в интернете и заменить их в датасете вручную."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['PRIM_LAT_DEC'] = df['PRIM_LAT_DEC'].replace(float('nan'), 38.892235)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['PRIM_LONG_DEC'] = df['PRIM_LONG_DEC'].replace(float('nan'), -77.003689)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Еще есть пропуски в целевой переменной Visitors. YearRaw у данных объектов тоже весьма странный - Total. Посмотрим на YearRaw."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['YearRaw'].value_counts().head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.loc[df['YearRaw'] == 'Total'].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "386 объектов имеют Total вместо года. Предполагается, что данные объекты были созданы в результате ошибки во время сбора данных. Предлагается удалить эти объекты."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.loc[df['YearRaw'] != 'Total']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Изменим формат целевой переменой Visitors и признака YearRaw."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['YearRaw'] = df['YearRaw'].astype('int64')\n",
    "df['Visitors'] = df['Visitors'].astype('int64')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Кроме того, для удобства переименуем YearRaw в Year."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.rename(columns={'YearRaw': 'Year'}, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Т.к. наша задача предсказать количество посетителей по историческим данных, большое количество исторических данных может затруднить задачу, в связи с чем облегчим себе задачу и оставим данные, начиная с 2000 года."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.loc[df['Year'] >= 2000]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### test и train"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь можно разбить данных не обучающую и тестовую выборки. Посмотрим, какой год последний в датасете."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sort(pd.unique(df['Year']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Будем предсказывать посетителей для 2016 года по данным с 2000 по 2015."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train = df.loc[df['Year'] != 2016]\n",
    "train.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test = df.loc[df['Year'] == 2016]\n",
    "test.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Обучающая выборка включает 5795 объектов, тестовая - 379. Всего 11 признаков."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Целевая переменная"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Сначала будем изучать целевую переменную Visitors. Проведем тесты на нормальность распределения и рассчитаем скошенность и эксцесс распределения."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(rc={\"figure.figsize\": (10, 6)})\n",
    "fig = plt.figure()\n",
    "sns.distplot(train['Visitors'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Из графика выше видно, что распределение совсем не похоже на нормальное. Чтобы убедиться в этом, добавим нормальное распределение (из scipy) к графику и сравним его с распределением целевой переменной."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(rc={\"figure.figsize\": (10, 6)})\n",
    "fig = plt.figure()\n",
    "sns.distplot(train['Visitors'], fit=stats.norm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Распределение целевой переменной скошенно и сильно отличается от нормального. Проведем стат.тесты."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Скошенность: %f\" % stats.skew(train['Visitors']))\n",
    "print(\"Эксцесс: %f\" % stats.kurtosis(train['Visitors']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.normaltest(df['Visitors'].values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как видим, скошенность и эксцесс распределения очень большие. В дальнейшем нам будет необходимо преобразовать целевую переменную и попытаться сделать распределение близким к нормальному."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Вещественные признаки"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Рассмотрим вещественный признаки."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В датасете всего три вещественных признака, построим графики, чтобы найти какие-нибудь закономерности или зависимости."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "visitors_df = train[['Visitors'] + ['Year']]\n",
    "visitors_df.groupby('Year').sum().plot(figsize=(12,6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Шкала графика выше вводит заблуждение, на самом деле количество посетителей меняется не так сильно от года в год. Чтобы убедиться в этом, посмотрим на bar-chart."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "visitors_df.groupby('Year').sum().plot(kind='bar', rot=0,figsize=(12,6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Что и требовалось доказать. Тем не менее, число посителей все же выросло в 2015 году. Теперь посмотрим на корреляции между веществеными признаками и целевой переменной."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(rc={\"figure.figsize\": (12, 6)})\n",
    "sns.heatmap(train.corr(), vmin=0, vmax=1,annot=True, cmap=\"Blues\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как видно из графика, корреляция между вещественными признаками и целевой переменной достаточно слабая."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Изучим датасет на выбросы. Для этого построим пару графиков и попытаемся найти и интерпретировать выбросы."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(rc={\"figure.figsize\": (12, 6)})\n",
    "sns.violinplot(x='Year', y=\"Visitors\", data=train, palette=\"muted\", split=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(rc={\"figure.figsize\": (12, 6)})\n",
    "plt.plot(train['Year'], train['Visitors'], 'r.')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Графики выше показывают, что у датасета есть выбросы, но необходимо изучить природу таких больших значений."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train.sort_values(by='Visitors', ascending=False).head(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Самое большое количество посещений на протяжении всего периода 2000-2015 у Parkway, что вполне логично. Таким образом, большие значения объясняются природой данных."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим на другие достопримечательности, которые собирают большое количество посетителей."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.unique(train['Unit Name'].loc[(train['Visitors'] >= 5000000) & (train['Unit Name'] != 'Blue Ridge Parkway')])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Названия объектов говорят сами за себя, они действительно набирают миллионы посетителей каждый год, поэтому причина больших значений вполне объяснима - природа данных."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Категориальные признаки"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В датасете 6 категориальных признаков. Один Parkname из них имеет пропуски, и кроме того, полные названия объектов представлены в признаке Unit Name. Таким образом, предлагается удалить этот признак. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train.drop(['Parkname'], axis=1, inplace=True)\n",
    "test.drop(['Parkname'], axis=1, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь у нас 5 категориальных признаком. Посмотрим на количество уникальных категорий в каждом из них."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(len(pd.unique(df['Region'])))\n",
    "print(len(pd.unique(df['State'])))\n",
    "print(len(pd.unique(df['Unit Code'])))\n",
    "print(len(pd.unique(df['Unit Name'])))\n",
    "print(len(pd.unique(df['Unit Type'])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Предлагается, каждый такой признак преобразовать с помощью One-Hot кодирования."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 4. Закономерности, \"инсайты\", особенности данных"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- В датасете есть пропуски - все они были решены\n",
    "- В датасете были обнаружены ошибки, допущенные в ходе сбора данных - все они были также решены\n",
    "- В датасете были выявлены выбросы, которые объясняются природой данных\n",
    "- Распределение целевой переменной не нормальное - преобразуем его в части Предобработка данных\n",
    "- Между целевой переменной и другими вещественными признаками была выявлена слабая корреляция. Таким образом, при построении модели будем делать акцент на категориальные признаки\n",
    "- В датасете 5 категориальных признаков, предалагется преобразовать их с помощью One-Hot Encoding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 5-6. Выбор метрики и модели"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В качестве метрики будем использовать r2_score (r_squared) - коэффициент детерминации. Под данной метрикой подразумевают долю дисперсии зависимой переменной, которая объясняется рассматриваемой моделью зависимости, то есть объясняющими переменными (признаками). \n",
    "Выбор обусловлен:\n",
    "- задачей регресии, т.к. данная метрика является универсальной мерой и идеально подходит для данной задачи\n",
    "- спецификой и особенностями данных"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В качестве модели будем использовать:\n",
    "- Ridge, т.к. одна из универсальных простых линейных моделей, необходимых для восстановления регрессии\n",
    "- CatBoostRegressor, т.к. данный датасет включает в себя большое количество категориальных признаков, с которыми CatBoost хорошо справляется"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 7-8. Предобработка данных и создание новых признаков"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Некоторые признаки мы уже предобработали ранее. Далее предлагается удалить ненужный признак GNIS ID, т.к. сам по себе он не несет никакой информации, необходимой для модели, но с помощью него мы уже достали широту и долготу."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train.drop(['Gnis Id'], axis=1, inplace=True)\n",
    "test.drop(['Gnis Id'], axis=1, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Целевой признак"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как было отмечено ранее, распределение целевой переменной не прошло тест на нормальность. Поэтому необходимо попытаться сделать его близким к нормальному. Рассмотрим два метотда: логарифмирование и метод [Бокса-Кокса](http://www.machinelearning.ru/wiki/index.php?title=%D0%9C%D0%B5%D1%82%D0%BE%D0%B4_%D0%91%D0%BE%D0%BA%D1%81%D0%B0-%D0%9A%D0%BE%D0%BA%D1%81%D0%B0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Сначала применим первый метод - логарифмирование."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(rc={\"figure.figsize\": (10, 6)})\n",
    "fig = plt.figure()\n",
    "sns.distplot(np.log1p((df['Visitors'].values)), fit=stats.norm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Скошенность: %f\" % stats.skew(np.log1p((df['Visitors'].values))))\n",
    "print(\"Эксцесс: %f\" % stats.kurtosis(np.log1p((df['Visitors'].values))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.normaltest(np.log1p((df['Visitors'].values)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как видно из результатов теста на нормальность и графика, нам удалось сделать распределение более нормальным, но все еще присутствует скошенность и эксцесс. Применим второй метод - меттод Бокса-Кокса."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# добавим единицу к посетителям, т.к. для данного метода значения должны быть строго больше нуля \n",
    "stats.boxcox(df['Visitors'].values+1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(rc={\"figure.figsize\": (10, 6)})\n",
    "fig = plt.figure()\n",
    "sns.distplot(stats.boxcox(train['Visitors'].values+1)[0], fit=stats.norm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.normaltest(stats.boxcox(train['Visitors'].values+1)[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Скошенность: %f\" % stats.skew(stats.boxcox(train['Visitors'].values+1)[0]))\n",
    "print(\"Эксцесс: %f\" % stats.kurtosis(stats.boxcox(train['Visitors'].values+1)[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В данном случае, метод Бокса-Кокса был более эффективен, поэтому преобразуем нашу целевую переменную с помощью этого метода. Также отделим целевую переменную от других признаков."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, y_train = train.drop('Visitors', axis=1).values, stats.boxcox(train['Visitors'].values+1)[0]\n",
    "X_test, y_test = test.drop('Visitors', axis=1).values, stats.boxcox(test['Visitors'].values+1)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  OHE для категориальных признаков "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Будем использовать готовые выборки (X_train, X_test) для CatBoostRegressor, т.к. он сам преобразует категориальные признаки. Для Ridge преобразуем категориальные признаки с помощью One-Hot Encoding. Напишем функцию для преобразования, а затем применим ее ко всем категоримальным признакам."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ohe_encoder(column,dftype):\n",
    "    enc = OneHotEncoder(sparse=False)\n",
    "    label_encoder = LabelEncoder()\n",
    "\n",
    "\n",
    "    values = np.array(df['{}'.format(column)].values)\n",
    "    values_labeled = label_encoder.fit_transform(values)\n",
    "\n",
    "    enc.fit(values_labeled.reshape(len(values_labeled), 1))\n",
    "\n",
    "    if dftype == 'train':\n",
    "        df_labeled = label_encoder.transform(train['{}'.format(column)])\n",
    "    else:\n",
    "        df_labeled = label_encoder.transform(test['{}'.format(column)])\n",
    "    \n",
    "    \n",
    "    df_encoded = enc.transform(df_labeled.reshape(len(df_labeled), 1))\n",
    "    \n",
    "    return df_encoded"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "train_region_ohe = ohe_encoder(column='Region',dftype = 'train')\n",
    "train_state_ohe = ohe_encoder(column='State',dftype = 'train')\n",
    "train_ucode_ohe = ohe_encoder(column='Unit Code',dftype = 'train')\n",
    "train_uname_ohe = ohe_encoder(column='Unit Name',dftype = 'train')\n",
    "train_utype_ohe = ohe_encoder(column='Unit Type',dftype = 'train')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "test_region_ohe = ohe_encoder(column='Region',dftype = 'test')\n",
    "test_state_ohe = ohe_encoder(column='State',dftype = 'test')\n",
    "test_ucode_ohe = ohe_encoder(column='Unit Code',dftype = 'test')\n",
    "test_uname_ohe = ohe_encoder(column='Unit Name',dftype = 'test')\n",
    "test_utype_ohe = ohe_encoder(column='Unit Type',dftype = 'test')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Соединим все полученные признаки в два датасета train_stacked и test_stacked."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_stacked = np.hstack([train_region_ohe,train_state_ohe,\n",
    "                           train_ucode_ohe,train_uname_ohe,\n",
    "                           train_utype_ohe,\n",
    "                           train[['Year', 'PRIM_LAT_DEC','PRIM_LONG_DEC']].values])\n",
    "test_stacked = np.hstack([test_region_ohe,test_state_ohe,\n",
    "                           test_ucode_ohe,test_uname_ohe,\n",
    "                           test_utype_ohe,\n",
    "                           test[['Year', 'PRIM_LAT_DEC','PRIM_LONG_DEC']].values])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 9. Кросс-валидация, подбор параметров"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ridge Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для Ridge укажем различные значения для основного параметра alpha - порога регуляризации. Чем выше значение, тем сильнее регуляризация.\n",
    "Далее обучим GridSearchCV для выявления лучших параметров."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = {'alpha' : [1,0.1,1e2,1e3,1e4]}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "linreg = Ridge()\n",
    "kf = KFold(n_splits=3, shuffle=True, random_state=17)\n",
    "linreg_cv = GridSearchCV(linreg, param_grid=params, scoring='r2', cv=kf)\n",
    "linreg_cv.fit(train_stacked, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Best parameters for Ridge: ', linreg_cv.best_params_)\n",
    "print('R2 score for Ridge: ', round(linreg_cv.best_score_, 5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "GridSearchCV показал лучшее значение порога регуляризации 0.1, его и будем использовать на тестовой выборке."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### CatBoostRegressor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для CatBoostRegressor явно укажем какие из признаков являются категориальными (features). Кроме того, зададим различные значения для глубины дерева (depths) и количества итераций (iterations)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "features = [0,1,2,3,4]\n",
    "depths = [7,8,9,10]\n",
    "iterations = [250,300]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Напишем собственную кастомную версию кросс-валидации для CatBoost и сравним резльтуты различных комбинаций параметров."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "counter = 1\n",
    "for x in depths:\n",
    "    for y in iterations:\n",
    "        print('Starting  '+ str(counter) + ' training: depth - ' + str(x) + ', iterations - ' + str(y))\n",
    "        kf = KFold(n_splits=3,shuffle=True,random_state=17)\n",
    "        results = 0\n",
    "        fold = 1\n",
    "        for train_index, test_index in kf.split(X_train):\n",
    "            cat_train, cat_test = X_train[train_index], X_train[test_index]\n",
    "            train_labels, test_labels = y_train[train_index], y_train[test_index]\n",
    "\n",
    "\n",
    "            cat = CatBoostRegressor(depth=x,iterations=y,random_seed=17,logging_level='Silent')\n",
    "            cat.fit(cat_train, train_labels, cat_features=features)\n",
    "            cat_pred = cat.predict(cat_test)\n",
    "            \n",
    "            results += r2_score(test_labels, cat_pred)\n",
    "            print(fold, 'fold score is: ', r2_score(test_labels, cat_pred))\n",
    "            fold += 1\n",
    "        \n",
    "        print('Average score for',str(counter),'training is:',results/3)\n",
    "        print('Done...' + '\\n')\n",
    "        counter += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Лучший резлуьтат (0.8170822248344393) на кросс-валидации показала модель с параметрами depth - 10, iterations - 300. Именно их будем использовать на тестовой выборке."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Часть 10. Построение кривых валидации и обучения"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Построим кривые валидации для Ridge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,6))\n",
    "param_range = [1,0.1,1e2,1e3,1e4]\n",
    "train_scores, test_scores = validation_curve(\n",
    "    Ridge(), train_stacked, y_train, param_name='alpha', param_range=param_range,\n",
    "    cv=3, scoring=\"r2\")\n",
    "train_scores_mean = np.mean(train_scores, axis=1)\n",
    "train_scores_std = np.std(train_scores, axis=1)\n",
    "test_scores_mean = np.mean(test_scores, axis=1)\n",
    "test_scores_std = np.std(test_scores, axis=1)\n",
    "\n",
    "plt.title(\"Validation Curve with Ridge\")\n",
    "plt.xlabel(\"alpha\")\n",
    "plt.ylabel(\"Score\")\n",
    "plt.ylim(0.0, 1.1)\n",
    "lw = 2\n",
    "plt.semilogx(param_range, train_scores_mean, label=\"Training score\",\n",
    "             color=\"darkorange\", lw=lw)\n",
    "plt.fill_between(param_range, train_scores_mean - train_scores_std,\n",
    "                 train_scores_mean + train_scores_std, alpha=0.2,\n",
    "                 color=\"darkorange\", lw=lw)\n",
    "plt.semilogx(param_range, test_scores_mean, label=\"Cross-validation score\",\n",
    "             color=\"navy\", lw=lw)\n",
    "plt.fill_between(param_range, test_scores_mean - test_scores_std,\n",
    "                 test_scores_mean + test_scores_std, alpha=0.2,\n",
    "                 color=\"navy\", lw=lw)\n",
    "plt.legend(loc=\"best\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как мы видим из графика, кривые валидации с увеличением порога регуляризации alpha немного расходятся (справа-налево), это означает, что с увеличением alpha может наблюдаться переобучение. Теперь построим кривые обучения для Ridge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,6))\n",
    "\n",
    "train_sizes, train_scores, test_scores = learning_curve(\n",
    "        Ridge(), train_stacked, y_train, cv=3)\n",
    "train_scores_mean = np.mean(train_scores, axis=1)\n",
    "train_scores_std = np.std(train_scores, axis=1)\n",
    "test_scores_mean = np.mean(test_scores, axis=1)\n",
    "test_scores_std = np.std(test_scores, axis=1)\n",
    "plt.grid()\n",
    "\n",
    "plt.fill_between(train_sizes, train_scores_mean - train_scores_std,\n",
    "                     train_scores_mean + train_scores_std, alpha=0.1,\n",
    "                     color=\"r\")\n",
    "plt.fill_between(train_sizes, test_scores_mean - test_scores_std,\n",
    "                     test_scores_mean + test_scores_std, alpha=0.1, color=\"g\")\n",
    "plt.plot(train_sizes, train_scores_mean, 'o-', color=\"r\",\n",
    "             label=\"Training score\")\n",
    "plt.plot(train_sizes, test_scores_mean, 'o-', color=\"g\",\n",
    "             label=\"Cross-validation score\")\n",
    "\n",
    "\n",
    "plt.title(\"Validation Curve with Ridge\")\n",
    "plt.xlabel(\"Training examples\")\n",
    "plt.ylabel(\"Score\")\n",
    "plt.legend(loc=\"best\")\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "На графике видно, что с увеличением количества примеров. качества на валидации улучшается. Т.к. кривые обучения не сблизились в самом конце, необходимо добавить больше примеров обуччающей выборки."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Часть 11. Прогноз для тестовой выборки и Оценка модели с описанием выбранной метрики"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь, когда у нас есть лучшие параметры, мы готовы обучить модели на тестовой выборке, т.е. предсказать количество посетителей национальных достопримечательностей США."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ridge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "ridge = Ridge(alpha=0.1)\n",
    "ridge.fit(train_stacked, y_train)\n",
    "ridge_pred = ridge.predict(test_stacked)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r2_score(y_test, ridge_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Коэффициент детерминации достаточно высокий но сильно отличается от результата, показанного на кросс-валидации (0.95302)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### CatBoostRegressor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "cat = CatBoostRegressor(depth=10, iterations=300,eval_metric='R2',logging_level='Silent')\n",
    "cat.fit(X_train, y_train, cat_features=features)\n",
    "cat_pred = cat.predict(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r2_score(y_test, cat_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CatBoostRegressor показал себя хуже на тестовой выборке (0.7529929633300172), чем лучшая модель на валидации (0.8170822248344393). Вполне вероятно, что виной такой разницы может быть переобучение."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Часть 12. Выводы "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В результате исследования, были построены две модели Ridge из sklearn и CatBoostRegressor из CatBoost.\n",
    "На тестовой выборке CatBoostRegressor (0.7529929633300172) показал результат хуже, чем Ridge (0.8711647097580936). Обе модели переобучаются. Если Ridge переобучается из-за высокого порога alpha, то для выяснения причин переобучения CatBoostRegressor потребуется более детальное исследование.\n",
    "Кривые обучения также показали, что необходимо добавить больше примеров для обучающей выборки.\n",
    "\n",
    "Предполагаемые улучшения:\n",
    "1. Избавление от переобучения\n",
    "2. Создание новых признаков из текущего датасета + можно привлечь сторонние дополнительные данные\n",
    "3. Тюнинг параметров, т.к. были расмотренны только основные параметры моделей.\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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}