# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# + [markdown] id="TK7oX5qdEEYr"
# <p style="align: center;"><img align=center src="https://s8.hostingkartinok.com/uploads/images/2018/08/308b49fcfbc619d629fe4604bceb67ac.jpg" width=500 height=450/></p>
#
# <h3 style="text-align: center;"><b>Школа глубокого обучения ФПМИ МФТИ</b></h3>
#
# <h3 style="text-align: center;"><b>Домашнее задание. Продвинутый поток. Весна 2021</b></h3>
#
# Это домашнее задание будет посвящено полноценному решению задачи машинного обучения.

# + [markdown] id="zlB-owfaEEYs"
# Есть две части этого домашнего задания: 
# * Сделать полноценный отчет о вашей работе: как вы обработали данные, какие модели попробовали и какие результаты получились (максимум 10 баллов). За каждую выполненную часть будет начислено определенное количество баллов.
# * Лучшее решение отправить в соревнование на [kaggle](https://www.kaggle.com/c/advanced-dls-spring-2021/) (максимум 5 баллов). За прохождение определенного порогов будут начисляться баллы.
#
#
# **Обе части будут проверяться в формате peer-review. Т.е. вашу посылку на степик будут проверять несколько других студентов и аггрегация их оценок будет выставлена. В то же время вам тоже нужно будет проверить несколько других учеников.**
#
# **Пожалуйста, делайте свою работу чистой и понятной, чтобы облегчить проверку. Если у вас будут проблемы с решением или хочется совета, то пишите в наш чат в телеграме или в лс @runfme. Если вы захотите проаппелировать оценку, то пипшите в лс @runfme.**
#
# **Во всех пунктах указания это минимальный набор вещей, которые стоит сделать. Если вы можете сделать какой-то шаг лучше или добавить что-то свое - дерзайте!**

# + [markdown] id="Tu_JvqcBEN8Y"
# # Как проверять?
#
# Ставьте полный балл, если выполнены все рекомендации или сделано что-то более интересное и сложное. За каждый отсустствующий пункт из рекомендация снижайте 1 балл.

# + [markdown] id="ninJ63mJEEYt"
# # Метрика
#
# Перед решением любой задачи важно понимать, как будет оцениваться ваше решение. В данном случае мы используем стандартную для задачи классификации метрику ROC-AUC. Ее можно вычислить используя только предсказанные вероятности и истинные классы без конкретного порога классификации + она раотает даже если классы в данных сильно несбалансированны (примеров одного класса в десятки раз больше примеров длугого). Именно поэтому она очень удобна для соревнований.
#
# Посчитать ее легко:
#

# + id="SQIrka7yEEYu" outputId="ef45a9f2-7571-47b1-e697-b474b65cb0fb"
from sklearn.metrics import roc_auc_score

y_true = [
    0,
    1,
    1,
    0,
    1
]

y_predictions = [
    0.1,
    0.9,
    0.4,
    0.6,
    0.61
]

roc_auc_score(y_true, y_predictions)

# + [markdown] id="IrDNNkNTEEYz"
# # Первая часть. Исследование

# + id="lzLqEeZKEEYz"
import pandas as pd
import numpy as np
import seaborn as sns

# %matplotlib inline
from matplotlib import pyplot as plt
plt.style.use('ggplot')

import warnings
warnings.simplefilter('ignore')

np.random.seed(42)

# + [markdown] id="LOlxdURSEEY3"
# ## Загрузка данных (2 балла)
#
# 1) Посмотрите на случайные строчки. 
#
# 2) Посмотрите, есть ли в датасете незаполненные значения (nan'ы) с помощью data.isna() или data.info() и, если нужно, замените их на что-то. Будет хорошо, если вы построите табличку с количеством nan в каждой колонке.

# + id="pw-Brue9EEY3"
data = pd.read_csv('./train.csv')
# -

# ### Краем глаза посмотрим на основные сводные таблицы по нашим данным (так называемый, первичный анализ данных...)

data.sample(n=10)  # слуайные объекты из тренировочной выборки

data.describe()  # основные статистики по числовым признакам

data.describe(include=('object'))  # оснонвые статистики по категориальным признакам

data.isna().sum() / data.shape[0] * 100  # смотрим на пропуски в процентах

# отрисуем график пропусков (синий - пропуск)
colours = ['#000099', '#ffff00'] 
sns.heatmap(
    data.isnull(),
    cmap=sns.color_palette(colours),
    yticklabels=False,
);

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

data.info()

# Как можно заметить, достаточно большое число данных у нас относится к категорильным, что, конечно, затруднит работу алгоритмов. Но ничего, это частично поправимо!

# + id="KgnkkF5bEEY9"
# Для вашего удобства списки с именами разных колонок

# Числовые признаки
num_cols = [
    'ClientPeriod',
    'MonthlySpending',
    'TotalSpent',
]

# Категориальные признаки
cat_cols = [
    'Sex',
    'IsSeniorCitizen',
    'HasPartner',
    'HasChild',
    'HasPhoneService',
    'HasMultiplePhoneNumbers',
    'HasInternetService',
    'HasOnlineSecurityService',
    'HasOnlineBackup',
    'HasDeviceProtection',
    'HasTechSupportAccess',
    'HasOnlineTV',
    'HasMovieSubscription',
    'HasContractPhone',
    'IsBillingPaperless',
    'PaymentMethod'
]

feature_cols = num_cols + cat_cols
target_col = 'Churn'

bin_cols = [  # добавим еще лист бинарных и тернарных признаков, но обзовем его бинарарным 
    'Sex',
    'IsSeniorCitizen',
    'HasPartner',
    'HasChild',
    'HasPhoneService',
    'HasMultiplePhoneNumbers',
    'HasInternetService',
    'HasOnlineSecurityService',
    'HasOnlineBackup',
    'HasDeviceProtection',
    'HasTechSupportAccess',
    'HasOnlineTV',
    'HasMovieSubscription',
    'HasContractPhone',
]

# + [markdown] id="dK370bPCEEZD"
# ## Анализ данных (3 балла)
#
# 1) Для численных призанков постройте гистограмму (*plt.hist(...)*) или boxplot (*plt.boxplot(...)*). Для категориальных посчитайте количество каждого значения для каждого признака. Для каждой колонки надо сделать *data.value_counts()* и построить bar диаграммы *plt.bar(...)* или круговые диаграммы *plt.pie(...)* (хорошо, елси вы сможете это сделать на одном гарфике с помощью *plt.subplots(...)*). 
#
# 2) Посмотрите на распределение целевой переменной и скажите, являются ли классы несбалансированными.
#
# 3) (Если будет желание) Поиграйте с разными библиотеками для визуализации - *sns*, *pandas_visual_analysis*, etc.
#
# Второй пункт очень важен, потому что существуют задачи классификации с несбалансированными классами. Например, это может значить, что в датасете намного больше примеров 0 класса. В таких случаях нужно 1) не использовать accuracy как метрику 2) использовать методы борьбы с imbalanced dataset (обычно если датасет сильно несбалансирован, т.е. класса 1 в 20 раз меньше класса 0).
# -

# ### Сперва узнаем, какое вообще "распределение" целевой переменной

# +
types, frequency = np.unique(data['Churn'], return_counts=True)
print(f"\t{round(frequency[0] / data['Churn'].shape[0], 2)} - rate of {types[0]}\n\
    \t{round(frequency[1] / data['Churn'].shape[0], 2)} - rate of {types[1]}")

data['Churn'].plot(kind='hist');


# -

# Видно, что выборка несбалансированная, так как объекты с целевой переменной `1` составляют лишь примерно четверть всех данных, но ресэмплить я все-таки не буду, попробую использовать другие методы просто более аккуратно. Например, не будем оценивать качество модели по `accuracy`
#
# Также стоит обратить внимание, что минимально допустимая точность на хорошем алгоритме должна быть минимум выше, чем `0.74`, так как такого качества возможно достичь даже "константным" алгоритмом с ответом `0`

# + [markdown] id="NZkbgFJZEEZE"
# ### Посмотрим теперь на распределение числовых признаков
# -

# #### Но сперва...

# Стоит обратить внимание, что почему-то `TotalSpent` имеет тип `object`, хотя и является числовым признаком. Также `IsSeniorCitizen` является чилсловым, хотя по сути это категориальный. Исправим же это недоразумение!
#
# В процессе кастования `TotalSpent` к `float64` выясняется, что  появились пропуски. Заполним их средним.

def prepare_data(dataset):
    new = dataset.copy()
    
    new['TotalSpent'] = pd.to_numeric(new['TotalSpent'], errors='coerce')
    new['IsSeniorCitizen'] = new['IsSeniorCitizen'].astype('object')
    new['TotalSpent'] = new['TotalSpent'].fillna(new['TotalSpent'].mean())
    
    assert data.isnull().sum().sum() == 0
    return new


data = prepare_data(data)

data.isnull().sum()

data.dtypes

# #### вернемся теперь к нашим числовым признакам

sns.pairplot(data[num_cols + [target_col]], hue='Churn');

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

data[num_cols].hist(figsize=(12, 8));

# То же самое, но в логарифмической шкале:

data[num_cols].hist(figsize=(12, 8), log=True);


# Видно, что в целом нельзя говорить о тяжелых хвостах в распределениях числовых признаков, но, тем не менее не побрезгуем затем проналазировать данные на предмет нахождения их в *окрестности трех $\sigma$* среднего

def detect_outlier_std(data, feature, threshold=3):
    mean = data[feature].mean()
    sigma = data[feature].std()
    
    upper_bound = mean + threshold * sigma
    lower_bound = mean - threshold * sigma
    
    outliers = pd.concat([data[data[feature] > upper_bound], data[data[feature] < lower_bound]], axis=0)
    
    return outliers


print("Согласно правилу трех сигм:\n")
for feature in num_cols:
    print(f"количество выбросов для признака {feature}: {detect_outlier_std(data, feature).shape[0]}")


# Как можно заметить, все числовые признаки лежат в окрестности трех сигм (трех стандартных отклонений). Аналогично, проделаем процедуру с интерквантильным размахом:

def detect_outlier_iqr(data, feature, threshold=1.5):
    IQR = data[feature].quantile(0.75) - data[feature].quantile(0.25)
    
    upper_bound = data[feature].quantile(0.75) + (IQR * threshold)
    lower_bound = data[feature].quantile(0.25) - (IQR * threshold)
    
    outliers = pd.concat([data[data[feature] > upper_bound], data[data[feature] < lower_bound]], axis=0)
    
    return outliers


print(f"Согласно правилу 1.5 интерквантильных размахов (1.5*IQR):\n")
for feature in num_cols:
    print(f"количество выбросов для признака {feature}: {detect_outlier_iqr(data, feature).shape[0]}")

# C чистой совестью заявляем, что явных выбросов в данных **НЕТ**

# Посмотрим еще на ящики с усами для числовых признаков относительно разных категориальных переменных

sns.boxplot(x='Sex', y='ClientPeriod', data=data, hue='HasPartner', saturation=1);

# Например, здесь видно, что если у человека есть партнер (вне зависимости от пола), то `ClientPeriod` больше, чем у людей без партнера

sns.boxplot(x='Churn', y='TotalSpent', hue='HasMovieSubscription',data=data);

# Из этого боксплота ничего особо интересного нет, так как он лишь говорит о том, что люди с подпиской тратят больше, чем люди без интернета или без подписки. Но пусть он все равно тут будет :)

# Также изобразим, как признаки связаны с оттоком

# +
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(16, 7))

for idx, feat in  enumerate(num_cols):
    sns.boxplot(x='Churn', y=feat, data=data, ax=axes[idx])
    axes[idx].set_xlabel('Churn')
    axes[idx].set_ylabel(feat);
# -

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

# Теперь поизучаем матрицу корреляций:

# +
fig, ax = plt.subplots(figsize=(7,6))

sns.heatmap(
    data[num_cols + [target_col]].corr(),
    annot=True,
);
# -

# Как видно из матрицы, `TotalSpent` и `ClientPeriod` довольно прилично коррелируют между собой. Тем не менее, убирать я не буду один из этих столбцов, так как коррелируют они не очень сильно. Заметим также, что `TotalSpent` и `MonthlySpending` имеют некоторую линейную зависимость, вероятно, это из-за того, что оба эти признака суть признаки про траты

# ### Перейдем к категориальным

# изобразим для бинарных (или почти) признаков `countplot` относительно целевой переменной

# +
fig, axes = plt.subplots(nrows=7, ncols=2, figsize=(16, 14))
plt.subplots_adjust(top=2)

for ax, binary_feature in zip(axes.flat, bin_cols):
    sns.countplot(x=binary_feature, hue='Churn', data=data, ax=ax)
    ax.set_xlabel('Churn')
    ax.set_ylabel(f'{binary_feature} count');
# -

# В целом видно, что среди тех, у кого нет интернет сервисов, довольно мало покидают.

# + [markdown] id="zg60u3QDEEZH"
# (Дополнительно) Если вы нашли какие-то ошибки в данных или выбросы, то можете их убрать. Тут можно поэксперементировать с обработкой данных как угодно, но не за баллы.

# + [markdown] id="fwfksF1gEEZI"
# А вроде бы их и нет...

# + [markdown] id="DviiJd8REEZK"
# ## Применение линейных моделей (3 балла)
#
# 1) Обработайте данные для того, чтобы к ним можно было применить LogisticRegression. Т.е. отнормируйте числовые признаки, а категориальные закодируйте с помощью one-hot-encoding'а. 
#
# 2) С помощью кроссвалидации или разделения на train/valid выборку протестируйте разные значения гиперпараметра C и выберите лучший (можно тестировать С=100, 10, 1, 0.1, 0.01, 0.001) по метрике ROC-AUC. 
#
# Если вы разделяете на train/valid, то используйте LogisticRegressionCV. Он сам при вызове .fit() подберет параметр С. (не забудьте передать scroing='roc_auc', чтобы при кроссвалидации сравнивались значения этой метрики, и refit=True, чтобы при потом модель обучилась на всем датасете с лучшим параметром C). 
#
#
# (более сложный вариант) Если вы будете использовать кроссвалидацию, то преобразования данных и LogisticRegression нужно соединить в один Pipeline с помощью make_pipeline, как это делалось во втором семинаре. Потом pipeline надо передать в GridSearchCV. Для one-hot-encoding'a можно испльзовать комбинацию LabelEncoder + OneHotEncoder (сначала превращаем строчки в числа, а потом числа првращаем в one-hot вектора.)

# + id="xHCLDmwqEEZL"
from sklearn.linear_model import LogisticRegression, LogisticRegressionCV
from sklearn.model_selection import train_test_split, GridSearchCV, RandomizedSearchCV
from sklearn.preprocessing import StandardScaler, RobustScaler, LabelEncoder, OneHotEncoder
from sklearn.pipeline import make_pipeline
from sklearn.metrics import roc_auc_score, make_scorer
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.compose import ColumnTransformer

from sklearn import set_config
set_config(display='diagram')

# + [markdown] id="P2Yv3uYtEEZO"
# ### Для начала разберемся с категориальными признаками
#
# Проведем операцию one-hot кодирования

# +
dummies = pd.get_dummies(data[cat_cols])

X = pd.concat([data[num_cols], dummies], axis=1)
y = data[target_col]
X.head()
# -

# ### Теперь отнормируем числовые признаки
#
# Сперва разобьем данные на train / test, чтобы не произошло вдруг внезапной утечки при нормировке по всему датасету

X_train, X_test, y_train, y_test = train_test_split(X.values, y.values, test_size=0.33)

# +
scaler = StandardScaler()

X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
# -

# ### Данные готовы к бою! \*Логистическая регрессия выходит в октагон...\*

# +
clf_vanilla = LogisticRegressionCV(
    Cs=(100, 10, 1, 0.1, 0.01, 0.001),
    cv=5,
    scoring='roc_auc',
    random_state=42,
    refit=True,
)

clf_vanilla.fit(X_train_scaled, y_train)

# +
y_train_prediction = clf_vanilla.predict_proba(X_train_scaled)[:, 1]

print(f"ROC-AUC on the train data: {roc_auc_score(y_train, y_train_prediction)}")

# +
y_test_prediction = clf_vanilla.predict_proba(X_test_scaled)[:, 1]

print(f"ROC-AUC on the test data: {roc_auc_score(y_test, y_test_prediction)}")
# -

print(f"the best coef of regulazation: {float(clf_vanilla.C_)}")

# ### Реализуем также более трудоемкий вариант

# +
num_transformer = Pipeline(
    steps=[
        ('imputer', SimpleImputer(strategy='mean')),
        ('scaler', StandardScaler()),
    ]
)

cat_tranformer = Pipeline(
    steps=[
        ('to_one_hot', OneHotEncoder()),
    ]
)

preprocessor = ColumnTransformer(
    transformers=[
        ('num', num_transformer, num_cols),
        ('cat', cat_tranformer, cat_cols)
    ]
)

clf = Pipeline(
    steps=[
        ('preprocessor', preprocessor),
        ('classifier', LogisticRegression())
    ]
)
clf
# -

# Заново разобьем сырые данные на train/test, чтобы убедиться что pipeline работает корректно

X_train, X_test, y_train, y_test = train_test_split(data[feature_cols], data[target_col], test_size=0.2)

# Теперь создаем `GridSearch`

# +
param_grid = {
    "classifier__C": (100, 10, 1, 0.1, 0.01, 0.001),
    "classifier__penalty": ["l1","l2"],
}

grid_search = GridSearchCV(clf, param_grid, cv=10)
grid_search.fit(X_train, y_train)
# -

print(("best logistic regression from grid search: %.3f"
       % grid_search.score(X_test, y_test)))

# + [markdown] id="BlTeVy7fEEZR"
# ## Применение градиентного бустинга (2 балла)
#
# Если вы хотите получить баллы за точный ответ, то стоит попробовать градиентный бустинг. Часто градиентный бустинг с дефолтными параметрами даст вам 80% результата за 0% усилий.
#
# Мы будем использовать catboost, поэтому нам не надо кодировать категориальные признаки. catboost сделает это сам (в .fit() надо передать cat_features=cat_cols). А численные признаки нормировать для моделей, основанных на деревьях не нужно.
#
# 1) Разделите выборку на train/valid. Протестируйте catboost cо стандартными параметрами.
#
# 2) Протестируйте разные занчения параметроа количества деревьев и learning_rate'а и выберите лучшую по метрике ROC-AUC комбинацию. 
#
# (Дополнительно) Есть некоторые сложности с тем, чтобы использовать CatBoostClassifier вместе с GridSearchCV, поэтому мы не просим использовать кроссвалидацию. Но можете попробовать)

# +
# # !pip install catboost
# -

from catboost import CatBoostClassifier, Pool

# +
X = data.drop(['Churn'], axis=1)
y = data['Churn']

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, shuffle=False)
cat_indx = data.columns.get_indexer(cat_cols)
# -

X_train.head()

# +
# grid_clf = CatBoostClassifier(verbose=False)

# params = {
#    "eval_metric": ['AUC'],
#    "use_best_model": [True],
#    "iterations": [ 800, 1000, 1200],
#   "learning_rate": [0.02, 0.04,],
#    "cat_features": [cat_indx],
#    "depth": [6, 8, 10],
#    "loss_function": ['Logloss', "CrossEntropy"],
#    'random_seed': [42],
#}

# catboost_grid = RandomizedSearchCV(estimator=grid_clf, param_distributions=params, cv=5)

# +
catboost_clf = CatBoostClassifier()

params = {
    'iterations': [1000],
    'depth': [6, 8, 10],
    'loss_function': ['Logloss', 'CrossEntropy'],
    'random_seed': [42],
    'verbose': [False],
    'cat_features': [cat_indx],
}

scorer = make_scorer(roc_auc_score, needs_proba=True)
catboost_grid = GridSearchCV(estimator=catboost_clf, param_grid=params, scoring=scorer, cv=5)
# -

catboost_grid.fit(X_train, y_train)

best_param = catboost_grid.best_params_
best_param

final_catboost_clf = CatBoostClassifier(**best_param,
                                        eval_metric='AUC',
                                        use_best_model=True,
                                        learning_rate=0.0499,
                                       )
train_pool =  Pool(X_train, y_train, cat_features=cat_indx)
final_catboost_clf.fit(train_pool, eval_set=(X_test, y_test));

# +
train_pred = final_catboost_clf.predict_proba(X_train)[:, 1]
test_pred = final_catboost_clf.predict_proba(X_test)[:, 1]


print(f"ROC-AUC on the train data: {roc_auc_score(y_train, train_pred)}")
print(f"ROC-AUC on the test data: {roc_auc_score(y_test, test_pred)}")

# + [markdown] id="dDMXbvNZEEZV"
# # Предсказания

# + id="G_P4wFNaEEZW" outputId="1fba5dfc-88e4-49e3-ed8a-afe21ae3325a"
best_model = final_catboost_clf
print(f"preferable number of tree : {best_model.tree_count_}")
# -

X_test = pd.read_csv('./test.csv')
X_test = prepare_data(X_test)

X_test.isnull().sum()

# + id="FfSufx0CEEZZ"
submission = pd.read_csv('./submission.csv')
predicted =  best_model.predict_proba(X_test)[:, 1]
submission['Churn'] =  predicted

submission.to_csv('./my_submission.csv', index=False)
# -

submission

# ## Полученный в итоге максимальный скор: 0.85279
# ## В соревновании ищите меня как `<NAME>`
#

# + [markdown] id="mkxjbGkVEEZc"
# # Kaggle (5 баллов)
#
# Как выставить баллы:
#
# 1) 1 >= roc auc > 0.84 это 5 баллов
#
# 2) 0.84 >= roc auc > 0.7 это 3 балла
#
# 3) 0.7 >= roc auc > 0.6 это 1 балл
#
# 4) 0.6 >= roc auc это 0 баллов
#
#
# Для выполнения задания необходимо выполнить следующие шаги.
# * Зарегистрироваться на платформе [kaggle.com](kaggle.com). Процесс выставления оценок будет проходить при подведении итогового рейтинга. Пожалуйста, укажите во вкладке Team -> Team name свои имя и фамилию в формате Имя_Фамилия (важно, чтобы имя и фамилия совпадали с данными на Stepik).
# * Обучить модель, получить файл с ответами в формате .csv и сдать его в конкурс. Пробуйте и экспериментируйте. Обратите внимание, что вы можете выполнять до 20 попыток сдачи на kaggle в день.
# * После окончания соревнования отправить в итоговый ноутбук с решением на степик. 
# * После дедлайна проверьте посылки других участников по критериям. Для этого надо зайти на степик, скачать их ноутбук и проверить скор в соревновании.
# -



StarcoderdataJupyterScriptsDedupFiltered