# ---
# 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
# ---

# # Домашняя работа 3. Задача классификации
#
# ~~*Дедлайн мягкий как облачко: 1 ноября, 21:00*~~
#
# ~~*Дедлайн жесткий как неспелая хурма: 5 ноября, 21:00*~~
#
# *Дедлайн унылый как Асино настроение: 6 ноября, 21:00* (жесткий и единственный, без снятия баллов)

# ### Оценивание и штрафы
#
# Максимальная оценка — 10 баллов. Еще есть 2 бонусных балла, которые можно добавить к любым домашкам или проверочным.
#
# Не списывайте, иначе всем участникам обнулим :)

# Для удобства проверки самостоятельно посчитайте свою максимальную оценку (исходя из набора решенных задач) и укажите ниже.
#
# **Оценка: 9.5**

print('Всем удачи!👒')

# +
from __future__ import annotations

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# %matplotlib inline

pd.set_option('display.max_rows', 100)
pd.set_option('display.max_columns', 100)
# -

# ## Часть 1. Логрег своими руками (4.5 балла)

# **Задание 1 (4 балла)**. Реализуйте логистическую регрессию, обучаемую с помощью:
# - градиентного спуска **(2 балла)**
#
# - стохастического градиентного спуска **(2 балла)**
#
# Во всех пунктах необходимо соблюдать два условия:
# - Циклы можно использовать только для итераций градиентного спуска;
# - В качестве критерия останова необходимо использовать (одновременно):
#
#     - проверку на евклидову норму разности весов на двух соседних итерациях (например, меньше некоторого малого числа порядка $10^{-6}$), задаваемого параметром `tolerance`;
#     - достижение максимального числа итераций (например, 10000), задаваемого параметром `max_iter`.
#
# Чтобы проследить, что оптимизационный процесс действительно сходится, добавим атрибут класса `loss_history`. В нём после вызова метода `fit` должны содержаться значения функции потерь для всех итераций градиентного спуска, начиная с нулевой.
#
# Инициализировать веса можно случайным образом или нулевым вектором.

# +
from sklearn.base import BaseEstimator

class LogReg(BaseEstimator):
    def __init__(self, gd_type: str = 'stochastic', tolerance: float = 1e-4, 
                 max_iter: int = 1000, eta: float = 1e-2,
                 w0: np.array = None) -> None:
        """
        Args:
          gd_type: Type of gradient descent ('full' or 'stochastic'). 
          tolerance: Threshold for stopping gradient descent.
          max_iter: Maximum number of steps in gradient descent.
          eta: Learning rate.
          w0: Array of shape d (d — number of weights to optimize).
              Initial weights.
        """
        self.gd_type = gd_type
        self.tolerance = tolerance
        self.max_iter = max_iter
        self.eta = eta
        self.w0 = w0
        self.w = None # was None
        self.loss_history = None 
    
    def fit(self, X: np.array, y: np.array) -> LogReg:
        """Fit the model on training data. Also, save value of loss after each iteration.
        Args:
          X: Training data.
          y: Target.
        Returns: 
          self: Fitted classsifier.  
        """
        X = np.c_[X, np.ones(y.size)]
        if self.w0 is None:
                  self.w0 = np.zeros(X.shape[1])
        
        self.w = self.w0

        if self.gd_type == "full":
            norm1 = np.linalg.norm(self.w)
            self.loss_history = []
            for _ in range(self.max_iter): # плюс self.calc или минус?
                self.w = self.w - self.calc_gradient(X, y)*self.eta   # /(1+i) # коэф 1/k
                norm2 = np.linalg.norm(self.w)
                self.loss_history = np.append(self.loss_history, self.calc_loss(X, y))
                if abs(norm2 - norm1) < self.tolerance:
                    break
                norm1 = norm2
            return self
        else:
            norm1 = np.linalg.norm(self.w)
            self.loss_history = []
            for _ in range(self.max_iter):
                i = np.random.randint(0, X.shape[0])
                Xi = X[i, :]
                yi = y[i]
                step = - self.calc_gradient(Xi, yi) * self.eta
                self.w += step
                norm2 = np.linalg.norm(self.w)
                if abs(norm2 - norm1) < self.tolerance:
                    return self
                norm1 = norm2
                self.loss_history = np.append(self.loss_history, self.calc_loss(X, y))
            return self
    
    
    def predict_proba(self, X: np.array) -> np.array:
        """Calculate probability of positive and negative class for each observation.
        Args:
          X: Array of shape (n, d).
             Data.
        Returns: 
             Array of shape (n, 2).
             Predicted probabilities.
        """
        X = np.c_[X, np.ones(X.shape[0])]
        if self.w is None:
            raise Exception('Not trained yet')
        return 1 / (1 + np.exp(-1 * (X.dot(self.w))))
        pass
    
    
    def predict(self, X: np.array) -> np.array:
        """Predict class for each observation.
        Args:
          X: Array of shape (n, d).
             Data.
        Returns: 
             Array of shape (n,).
             Predicted class labels.  
        """
        X = np.c_[X, np.ones(X.shape[0])]
        if self.w is None:
            raise Exception('Not trained yet')
        return np.sign(X.dot(self.w)) # может вернуть 0
        pass
    
    def calc_gradient(self, X: np.array, y: np.array) -> np.array:
        """Calculate gradient of loss function after each iteration.
        Args:
          X: Array of shape (n, d), n can be equal to 1 if 'stochastic'.
          y: Array of shape (n,).
        Returns:
          Array of shape (d,).
          Gradient of loss function after current iteration.
        """
        if self.gd_type == "stochastic":
            X = X.reshape(1, -1)
        M = y * (X.dot(self.w))
        a = y * np.exp(-M) / (1 + np.exp(-M))
        grad = -(X.T@ a) / y.size
        return grad


    def calc_loss(self, X: np.array, y: np.array) -> float:
        """Calculate value of loss function after each iteration.
        
        Args:
          X: Array of shape (n, d).
          y: Array of shape (n,).
        
        Returns:
          Value of loss function after current iteration.
        """ 
        return np.mean(np.log(1 + np.exp(-1 * X.dot(self.w))))
        pass
# -

# Далее предполагается, что вы используете собственную реализацию логистической регрессии.
# Если с написанием класса возникли проблемы, используйте реализацию sklearn, чтобы не терять баллы за остальные задания.
#
# В части 2 и далее я бы всем советовала использовать реализацию sklearn.

# Сгенерируем синтетические данные.

# +
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split

X, y = make_classification(
    n_samples=10000, n_features=10, n_informative=5, n_redundant=5,
    random_state=42)

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=42)
# -

# **Задание 2 (0.5 балла).** Обучите логистическую регрессию на синтетических данных. Нарисуйте изменение лосса во время обучения.



# +
cls1 = LogReg(max_iter=10000, tolerance=-1e-8)
cls1.fit(X_train, y_train)
cls1.predict(X_test)

cls2 = LogReg(gd_type="full",tolerance=1e-8, max_iter = 10000)
cls2.fit(X_train, y_train)
cls2.predict(X_test)

plt.plot(cls1.loss_history, label="stochastic")
plt.plot(cls2.loss_history, label="full")
plt.grid()
plt.legend()
# -

# На тестовой части посчитайте ROC-AUC, PR-AUC. Постройте ROC и PR кривые.

# +
#ROC

from sklearn import metrics
from sklearn.metrics import roc_curve
from sklearn.metrics import roc_auc_score
from sklearn.metrics import precision_recall_curve
from sklearn.metrics import average_precision_score

cls1 = LogReg(gd_type="full",tolerance=1e-8, max_iter = 10000)
cls1.fit(X_train, y_train)

cls2 = LogReg(max_iter=10000, tolerance=-1e-8)
cls2.fit(X_train, y_train)

y1_predicted =  cls1.predict_proba(X_test)
y2_predicted = cls2.predict_proba(X_test)

fig, ax = plt.subplots()
fpr1, tpr1, thr1 = roc_curve(y_test, y1_predicted)
fpr2, tpr2, thr2 = roc_curve(y_test, y2_predicted)

p1, r1, thresholds = precision_recall_curve(y_test, y1_predicted)
p2, r2, thresholds = precision_recall_curve(y_test, y2_predicted)

plt.plot(fpr1, tpr1, label="full")
plt.plot(fpr2, tpr2, label="stochastic")
plt.legend()
plt.title("ROC curves")
ax.set_xlabel('FPR')
ax.set_ylabel('TPR')
plt.grid()

# -

fig, ax = plt.subplots()
plt.plot(r1, p1, label="full")
plt.plot(r2, p2, label="stochastic")
plt.legend()
plt.title("PR curves")
ax.set_ylabel('Precision')
ax.set_xlabel('Recall')
plt.grid()

print("ROC AUC score (full):", roc_auc_score(y_test, y1_predicted))
print("ROC AUC score (stochastic):", roc_auc_score(y_test, y2_predicted))

ans1 = 0
ans2 = 0
ans1 = average_precision_score(y_test, y1_predicted)
ans2 = average_precision_score(y_test, y2_predicted)
print("PR AUC (full):", ans1)
print("PR AUC (stochastic):", ans2)

# ## Часть 2. Работа с категориальными признаками (2.5 балла)

# В этой части мы научимся обрабатывать категориальные переменные. Как вы уже знаете, закодировать их в виде столбика чисел недостаточно (это задаёт некоторый порядок, которого на категориальных переменных может и не быть, но модель попробует его выучить). Существует два основных способа обработки категориальных значений:
# - One-hot-кодирование
# - Счётчики (CTR, mean-target кодирование, ...) — каждый категориальный признак заменяется на среднее значение целевой переменной по всем объектам, имеющим одинаковое значение в этом признаке.
#
# Начнём с one-hot-кодирования. Допустим наш категориальный признак $f_j(x)$ принимает значения из множества $C=\{c_1, \dots, c_m\}$. Заменим его на $m$ бинарных признаков $b_1(x), \dots, b_m(x)$, каждый из которых является индикатором одного из возможных категориальных значений:
# $$
# b_i(x) = [f_j(x) = c_i]
# $$

# __Подготовка данных.__
#
# Загрузим данные [UCI Bank Marketing Dataset](https://archive.ics.uci.edu/ml/datasets/bank+marketing). Этот датасет содержит информацию о маркетинговой кампании какого-то банка, объектом в нем является телефонный звонок потенциальному клиенту с предложением некоторой услуги (утверждается, что это краткосрочный депозит), целевой переменной — ответ клиента (согласился ли он открыть депозит?). В качестве признакового описания используются характеристики клиента (образование, брак и т.д.), данные о звонке и различные экономические индикаторы — более подробная информация на страничке с датасетом.

# # !wget https://archive.ics.uci.edu/ml/machine-learning-databases/00222/bank-additional.zip
# # !unzip bank-additional.zip
df = pd.read_csv('C:/Users/vesel/Downloads/bank-additional-full.csv', sep=';')

y = df["y"]
y = (y == "yes").astype(np.int8)
df = df.drop(["y", "duration"], axis=1)

df.head(5)

df.describe()

# __Задание 3 (0.5 балла).__ Разделите выборку на обучающую и тестовую в соотношении 3:1. Зафиксируйте `random_state=777`, укажите значение параметра `stratify`. Один из столбцов (помимо таргета :) ) стоит сразу выкинуть из обучающей выборки. Какой? Не отказывайте себе. 
#
# Ответ: Нужно выбросить duration, этот показатель становится известен уже после совершения звонка, поэтому он не годится для предсказаний.
#

# df.groupby(["default", "campaign"]).count()
X_train, X_test, y_train, y_test = train_test_split(
    df, y, test_size=0.25, random_state=777, stratify=y)




# Закодируйте категориальные признаки с помощью `OrdinalEncoder`. Посчитайте качество (в этом задании будем работать c `AUC-PR`) при применении логистической регрессии. Здесь и далее для реализации последовательности этих действий (обработка признаков + обучение модели) используйте пайплайны. Замерьте время, потребовавшееся на обучение модели (с учетом кодирования признаков).
#
# __Вопрос__: почему в данном задании мы выбрали метрикой именно `AUC-PR`, а не, к примеру, `AUC-ROC`?
#
# __Ваш ответ__: Потому что в данном случае мы очень хотим угадывать тех, у кого высокая вероятность заключить с нами контракт

# +
numeric_features = df.select_dtypes([np.number]).columns
ctgs = df.columns.drop(numeric_features)

from sklearn.preprocessing import OrdinalEncoder
from sklearn.pipeline import make_pipeline
from sklearn.pipeline import Pipeline
from sklearn.linear_model import LogisticRegression
from sklearn.compose import ColumnTransformer, make_column_selector
from sklearn.preprocessing import StandardScaler


prep_ord = ColumnTransformer([
    ('enc', OrdinalEncoder(), ctgs),
    ('scaling', StandardScaler(), numeric_features)])

pipe_ord = Pipeline(steps=[('prep_ord' , prep_ord), ('clf', LogisticRegression(max_iter=2000))])
# enc = OrdinalEncoder()
# enc.fit(df[ctgs])
# enc.categories_
# df[ctgs] = enc.transform(df[ctgs])
# -

# %%timeit
pipe_ord.fit(X_train, y_train)
y_predicted = pipe_ord.predict_proba(X_test)
# print(y_predicted)




# +
def plotPR(y_predicted, lbl=""):
    p, r, thresholds = precision_recall_curve(y_test, y_predicted[:, 1])
    plt.plot(r, p, label=lbl)
    plt.title("PR curves")
    plt.ylabel('Precision')
    plt.xlabel('Recall')
    plt.grid()
    plt.legend()
    print(lbl ,"AUC-RP:", average_precision_score(y_test, y_predicted[:, 1]))

pipe_ord.fit(X_train, y_train)
y_predicted_ord = pipe_ord.predict_proba(X_test)

plotPR(y_predicted_ord, lbl="OrdinalEncoding")
# -







# __Задание 4 (0.5 балла).__ Закодируйте все категориальные признаки с помощью one-hot-кодирования. Обучите логистическую регрессию и посмотрите, как изменилось качество модели. Измерьте время, потребовавшееся на кодирование категориальных признаков и обучение модели.
#
#

# +
from sklearn.preprocessing import OneHotEncoder

X_train, X_test, y_train, y_test = train_test_split(
    df, y, test_size=0.25, random_state=777, stratify=y)

prep_onehot = ColumnTransformer([
    ('onehot', OneHotEncoder(), ctgs),
    ('scaling', StandardScaler(), numeric_features)])

pipe_onehot = Pipeline(steps=[('prep_onehot' , prep_onehot), ('clf', LogisticRegression(max_iter=2000))])

# -

# %%timeit
pipe_onehot.fit(X_train, y_train)
y_predicted = pipe_onehot.predict_proba(X_test)

# 3.11 секунд по сравнению с 0.735 секунд для 2000 итераций GD

pipe_onehot.fit(X_train, y_train)
y_predicted_onehot = pipe_onehot.predict_proba(X_test)

plotPR(y_predicted_onehot, lbl="OneHot")
plotPR(y_predicted_ord, lbl = "Ordinal")
# немного выросло

# +
#╰( ͡° ͜ʖ ͡° )つ──☆*:・ﾟ
# -

# Как можно было заменить, one-hot-кодирование может сильно увеличивать количество признаков в датасете, что сказывается на памяти, особенно, если некоторый признак имеет большое количество значений. Эту проблему решает другой способ кодирования категориальных признаков — счётчики. Основная идея в том, что нам важны не сами категории, а значения целевой переменной, которые имеют объекты этой категории. Каждый категориальный признак мы заменим средним значением целевой переменной по всем объектам этой же категории:
# $$
# g_j(x, X) = \frac{\sum_{i=1}^{l} [f_j(x) = f_j(x_i)][y_i = +1]}{\sum_{i=1}^{l} [f_j(x) = f_j(x_i)]}
# $$
#
# __Задание 5 (0.5 балла).__ Закодируйте категориальные переменные с помощью счётчиков (ровно так, как описано выше без каких-либо хитростей). Обучите логистическую регрессию и посмотрите на качество модели на тестовом множестве. Сравните время обучения с предыдущими экспериментами (с учетом кодирования признаков). Заметили ли вы что-то интересное?

X_train, X_test, y_train, y_test = train_test_split(
    df, y, test_size=0.25, random_state=777, stratify=y)

df = pd.read_csv('C:/Users/vesel/Downloads/bank-additional-full.csv', sep=';')
df["y"] = (df["y"] == "yes").astype(np.int8)


# +
# %%timeit

data=df.copy()
for i in ctgs:
        data = data.replace({f'{i}': dict(df.groupby(f'{i}')['y'].mean())})
data = data.drop(['y','duration'], axis=1)

X_train, X_test, y_train, y_test = train_test_split(
    data, y, test_size=0.25, random_state=777, stratify=y) # почему то без первых двух строк выдаёт ошибку


lg = LogisticRegression(max_iter=2000)
lg.fit(X_train, y_train)
y_predicted = lg.predict_proba(X_test)

# +

data=df.copy()
for i in ctgs:
        data = data.replace({f'{i}': dict(df.groupby(f'{i}')['y'].mean())})
data = data.drop(['y','duration'], axis=1)

X_train, X_test, y_train, y_test = train_test_split(
    data, y, test_size=0.25, random_state=777, stratify=y) # почему то без первых двух строк выдаёт ошибку
print(X_train)
lg = LogisticRegression(max_iter=2000)
lg.fit(X_train, y_train)
y_predicted_ctr = lg.predict_proba(X_test)
# -

plotPR(y_predicted_ctr)

plt.figure(figsize=(9, 9))
plotPR(y_predicted_ctr, lbl="Счётчики")
plotPR(y_predicted_ord, lbl="Ordinal")
plotPR(y_predicted_onehot, lbl="OneHot")

# Отметим, что такие признаки сами по себе являются классификаторами и, обучаясь на них, мы допускаем «утечку» целевой переменной в признаки. Это ведёт к переобучению, поэтому считать такие признаки необходимо таким образом, чтобы при вычислении счетчика для конкретного объекта его целевая метка не использовалась. Это можно делать следующими способами:
# 1. Вычислять значение счётчика по всем объектам расположенным выше в датасете (например, если у нас выборка отсортирована по времени).
# 2. Вычислять по фолдам, то есть делить выборку на некоторое количество частей и подсчитывать значение признаков по всем фолдам кроме текущего (как делается в кросс-валидации).
# 3. Вносить шум в посчитанные признаки. 
#
# __Задание 6 (0.5 балла).__ Реализуйте корректное вычисление счётчиков самым простым способом — добавлением шума к значениям (постарайтесь найти баланс между борьбой с переобучением и сохранением полезности признаков). Снова обучите логистическую регрессию, оцените качество. Сделайте выводы.
# Какие плюсы и минусы использования счётчиков по сравнению с one-hot-кодированием можно отметить?
#
# Ответ: С счётчиками df занимает меньше памяти и быстрее обучается, но обладает меньшим качеством

# +
#╰( ͡° ͜ʖ ͡° )つ──☆*:・ﾟ
X_train, X_test, y_train, y_test = train_test_split(
    data, y, test_size=0.25, random_state=777, stratify=y)


i_max = 0.01
score_max = 0
print(X_train)
for i in np.arange(0.001, 0.011, 0.001):
    X_train_с = X_train
    X_train_с[ctgs] += np.random.normal(0, i, size = X_train[ctgs].shape)
    
    scl = StandardScaler()
    X_train_с = scl.fit_transform(X_train_с.values)
    X_test = scl.transform(X_test)
    
    lg = LogisticRegression(max_iter=2000)
    lg.fit(X_train_с, y_train)
    y_predicted_noise = lg.predict_proba(X_test)
    score = average_precision_score(y_test, y_predicted_noise[:, 1])
    print(score)
    if score > score_max:
        i_max = i
        score_max = score
i_max

# +
i_max = 0.001
X_train, X_test, y_train, y_test = train_test_split(
    data, y, test_size=0.25, random_state=777, stratify=y)


X_train_с = X_train
X_train_с[ctgs] += np.random.normal(scale=i_max, size = X_train[ctgs].shape)

scl = StandardScaler()
X_train_с = scl.fit_transform(X_train_с)
X_test = scl.transform(X_test)

lg = LogisticRegression(max_iter=2000)
lg.fit(X_train_с, y_train)

y_predicted_noise = lg.predict_proba(X_test)
# -

plt.figure(figsize=(9, 9))
plotPR(y_predicted_ctr, lbl="Счётчики")
plotPR(y_predicted_noise, lbl="Счётчики с шумом")
plotPR(y_predicted_ord, lbl="Ordinal")
plotPR(y_predicted_onehot, lbl="OneHot")

# __Задание 7 (0.5 балла).__ В данных имеется признак «возраст клиента». Сейчас мы интерпретируем его как числовой, что в общем случае для линейной модели может быть неверной гипотезой. Тем не менее, у этого признака есть довольно много уникальных значений (сколько?), поэтому применять к нему one-hot кодирование может оказаться излишним. Попробуйте закодировать возраст с помощью счетчиков. Стало ли лучше?
#
#

# +

X_train, X_test, y_train, y_test = train_test_split(
    data, y, test_size=0.25, random_state=777, stratify=y)

ctgs_age = pd.Index(np.append( ["age"], ctgs))
data=df.copy()
for i in ctgs_age:
        data = data.replace({f'{i}': dict(df.groupby(f'{i}')['y'].mean())})
data = data.drop(['y','duration'], axis=1)

print(data)
X_train, X_test, y_train, y_test = train_test_split(
    data, y, test_size=0.25, random_state=777, stratify=y)

#╰( ͡° ͜ʖ ͡° )つ──☆*:・ﾟ
lg = LogisticRegression(max_iter=2000)
lg.fit(X_train, y_train)
y_predicted_age = lg.predict_proba(X_test)
# -

plt.figure(figsize=(9, 9))
plotPR(y_predicted_age, lbl="age in ctgs")
plotPR(y_predicted_ctr, lbl="Счётчики")
plotPR(y_predicted_noise, lbl="Счётчики с шумом")
plotPR(y_predicted_ord, lbl="Ordinal")
plotPR(y_predicted_onehot, lbl="OneHot") # СТАЛО СИЛЬНО ХУЖЕ

# Можно пойти и в обратную сторону. У нас есть признаки «месяц и день недели» для звонка. Попробуйте интерпретировать их как числовые (месяц от 0 до 12, дни недели от 0 до 4). Стало ли лучше в этот раз?

# +
#╰( ͡° ͜ʖ ͡° )つ──☆*:・ﾟ

data=df.copy()
for i in ctgs.drop(["month", "day_of_week"]):
        data = data.replace({f'{i}': dict(df.groupby(f'{i}')['y'].mean())})
data = data.drop(['y','duration'], axis=1)

X_train, X_test, y_train, y_test = train_test_split(
    df, y, test_size=0.25, random_state=777, stratify=y)
    

prep = ColumnTransformer([
    ('enc', OrdinalEncoder(), ["month", "day_of_week"]) ])

pipe = Pipeline(steps=[('prep' , prep), ('clf', LogisticRegression(max_iter=2000))])

pipe.fit(X_train, y_train)
y_predicted = pipe.predict_proba(X_test)

plotPR(y_predicted)
# -

# ## Часть 3. Отбор признаков (1 балл + 1 бонусный балл)

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

# +
df = pd.read_csv('C:/Users/vesel/Downloads/bank-additional-full.csv', sep=';')

X = df.drop(columns=['duration', 'y'])
y = (df.y == 'yes')


X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=777, stratify=y)
# -

# Как вы помните, в данных много категориальных признаков (сейчас давайте интерпретировать возраст как числовой). Давайте закодируем их с помощью one-hot кодирования. Исходные колонки с категориальными признаками можно удалить. Сколько признаков мы получили?



#╰( ͡° ͜ʖ ͡° )つ──☆*:・ﾟ
encoder = OneHotEncoder(sparse=False)
encoder.fit(X[ctgs])
a = encoder.transform(X[ctgs]) 
data = np.hstack((a, np.array(X.drop(ctgs, axis=1))))
data = pd.DataFrame(data)
# Получили 62 признака
data.shape



# Как и ранее, в качестве основной модели будем использовать логистическую регрессию, а целевой метрикой выберем `AUC-PR`. Обучите модель и посчитайте качество на тестовой выборке. Давайте запомним полученное значение.

# +
#╰( ͡° ͜ʖ ͡° )つ──☆*:・ﾟ
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=777, stratify=y)

categorical = list(X_train.dtypes[X_train.dtypes == "object"].index)
numeric_features = X_train.select_dtypes([np.number]).columns

column_transformer = ColumnTransformer([
    ('ohe', OneHotEncoder(), categorical),
    ('scaling', StandardScaler(), numeric_features)
])

pipeline = Pipeline(steps=[ ('ohe_and_scaling', column_transformer),
                           ('regression', LogisticRegression(penalty="none" ,max_iter=2000)) ])

model = pipeline.fit(X_train, y_train)
y_predicted = model.predict_proba(X_test)
    
plotPR(y_predicted)

# +

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=777, stratify=y)

categorical = list(X_train.dtypes[X_train.dtypes == "object"].index)
numeric_features = X_train.select_dtypes([np.number]).columns

column_transformer = ColumnTransformer([
    ('ohe', OneHotEncoder(), categorical),
])

pipeline = Pipeline(steps=[ ('ohe_and_scaling', column_transformer),
                           ('regression', LogisticRegression(penalty="none" ,max_iter=2000)) ])

model = pipeline.fit(X_train, y_train)
y_predicted = model.predict_proba(X_test)
    
plotPR(y_predicted)
# -

# Замечаем, что если делать без Scaling, то качество падает.

# ### Встроенные методы

# Допустим, мы хотим оставить только 40 лучших признаков. Попробуем сделать это несколькими способами.
#
# Начнём с отбора признаков с помощью линейной модели. Как известно, веса линейной модели можно интерпретировать как вклад каждого признака в предсказание таргета, а значит, модуль этого вклада можно интерпретировать как важность признаков. Такой метод отбора называются встроенным (embedded methods), так как он заложен в особенности модели.
#
# __Задание 8 (0.5 балла).__ Оставьте 40 признаков с наибольшим модулем соответствующего параметра линейной модели. Обучите модель заново и оцените её качество. Замерьте скорость такого отбора признаков.
#
# Изменилось ли качество? Как?
#
#

# +
#╰( ͡° ͜ʖ ͡° )つ──☆*:・ﾟ
encoder = OneHotEncoder(sparse=False)
encoder.fit(X[ctgs])
a = encoder.transform(X[ctgs]) 
data = np.hstack((a, np.array(X.drop(ctgs, axis=1))))
data = pd.DataFrame(data)

X_train, X_test, y_train, y_test = train_test_split(data, y, test_size=0.2, random_state=777, stratify=y)

model = LogisticRegression()
model.fit(X_train, y_train)
w = model.coef_[0]

edge = sorted(w, reverse = True, key=lambda x: abs(x))[39]
X_train = X_train.drop(X_train.columns[abs(w) < edge ], axis=1) # проверено, работает
X_test = X_test.drop(X_test.columns[abs(w) < edge ], axis=1) 


# y_predicted = model.predict_proba(X_test)
# plotPR(y_predicted)
model.fit(X_train, y_train)
y_pred = model.predict_proba(X_test)

plotPR(y_pred)
# -

# А теперь давайте подумаем, что мы не учли. Мы предположили, что признаки вносят вклад равномерно, но не учли их масштаба. Если мы умножим один из признаков в 100 раз, то без учёта регуляризации его вес уменьшится в эти же 100 раз. А мы на основе этого отбираем признаки! Давайте сначала отмасштабируем признаки каким-то известным вам способом, а только потом будем удалять признаки. 
#
# Кстати, в таком случае надо пересчитать качество на всех признаках (это тоже сделайте ниже).
#
# Что получилось?

# +
#╰( ͡° ͜ʖ ͡° )つ──☆*:・ﾟ
# Качество на всех признаках было подсчитано выше

X_train, X_test, y_train, y_test = train_test_split(data, y, test_size=0.2, random_state=777, stratify=y)
print(type(X_train))
scl = StandardScaler()
X_train = pd.DataFrame(scl.fit_transform(X_train))
X_test = pd.DataFrame(scl.transform(X_test))
print(type(X_train))

model = LogisticRegression()
model.fit(X_train, y_train)
w = model.coef_[0]

edge = sorted(w, reverse = True, key=lambda x: abs(x))[39]
X_train = X_train.drop(X_train.columns[abs(w) < edge ], axis=1) # проверено, работает
X_test = X_test.drop(X_test.columns[abs(w) < edge ], axis=1) 


# y_predicted = model.predict_proba(X_test)
# plotPR(y_predicted)
model.fit(X_train, y_train)
y_pred = model.predict_proba(X_test)
print(X_train.shape)
plotPR(y_pred)
# -

# ### Методы фильтрации
#
#
# Давайте отбирать признаки умнее, а именно через подсчёт некоторой функции для каждого признака. На основании значений этой функции будем оставлять наиболее важные признаки. Методы этого семейства называют фильтрующими или  filter methods. 
#
# В качестве такой функции будем считать t-статистику:
#
# $$t(j) = \frac{|\mu_+ - \mu_-|}{\sqrt{\frac{n_+ \sigma^2_+ + n_- \sigma^2_-}{n_+ + n_-}}},$$
#
# где $\mu$, $\sigma$, $n$ соответственно среднее, стандартное отклонение и количество объектов каждого из классов.
#
# __Задание 9 (0.5 балла).__ Оставьте 40 признаков с наибольшим значением $t$ и замерьте качество. Не забудьте замерить скорость отбора признаков в этом случае.
#
#

# +
#╰( ͡° ͜ʖ ͡° )つ──☆*:・ﾟ
# -

# ### Методы-обёртки
#
# __Бонус (1 балл).__ 
#
# Третий из рассматриваемых нами методов работает следующим образом: мы исключаем по очереди один из признаков и смотрим, как это влияет на качество. Удаляем признаки таким жадным способом, пока не окажется выполненым некоторое условие (количество признаков или ухудшение качества).
#
# Заметим, что во время отбора признаков нельзя подсматривать в тестовую выборку (так же как и при настройке гиперпараметров). Разделите обучающую выборку на 2 части, на одной из них обучайте модель без одного из признаков,  на второй части оценивайте качество. Исходную тестовую выборку стоит использовать только для финальной оценки качества.
#
# Снова оставьте только 40 признаков и оцените качество на тестовой выборке. Сколько времени занял такой отбор признаков?

# +
#╰( ͡° ͜ʖ ͡° )つ──☆*:・ﾟ
df = pd.read_csv('C:/Users/vesel/Downloads/bank-additional-full.csv', sep=';')

X = df.drop(columns=['duration', 'y'])
y = (df.y == 'yes')

encoder = OneHotEncoder(sparse=False)
encoder.fit(X[ctgs])
a = encoder.transform(X[ctgs]) 
data = np.hstack((a, np.array(X.drop(ctgs, axis=1))))
data = pd.DataFrame(data)

X_train, X_test, y_train, y_test = train_test_split(data, y, test_size=0.2, random_state=777, stratify=y)

# -

# Стоит отметить, что с помощью такого метода можно пойти и в обратную сторону. Попробуйте _добавлять_ по одному самому полезному признаку в выборку до тех пор, пока не наберется 40 штук. Найдется ли порог, при котором добавление следующих признаков будет только ухудшать качество модели?

X_train

#╰( ͡° ͜ʖ ͡° )つ──☆*:・ﾟ
cnt = 0
features = []
mx = 0
for i in range(62):
    features.append(i)
    model = LogisticRegression(max_iter=2000)
    model.fit(X_train[features], y_train)
    y_pred = model.predict_proba(X_test[features])
    score = average_precision_score(y_test, y_pred[:, 1])
#     print(y_pred)
    if score > mx:
        mx = score
    else:
        features.pop()
#     if len(features) == 40:
#         break
print(i, len(features), features)

model = LogisticRegression(max_iter=2000)
model.fit(X_train[features], y_train) # Можно оставить только 40 features[:40] но это снизит качество
y_pred = model.predict_proba(X_test[features])
plotPR(y_pred)

# Заметим, что алгоритм выбрал 44 признака, то есть он остановится только после выбора 40. Добавлять признаки можно было по другому, от способа добавления зависит какие признаки туда попадут.

# **Вопрос для всех**: Давайте подведём итоги по отбору признаков. Назовите преимущества и недостатки каждого из методов. Какой метод привёл к наилучшему качеству?
#
# **Ответ:** В основном onehot показывает лучшее качество, но может быть долгим и занимать много памяти, Счётчики также работают неплохо но склонны переобучаться, при поправках, например, шумах работают неплохо, также они работают достаточно быстро - быстрее onehot. Остальные проигрывают по первым двум, хоть и могут быть быстрее.

# ## Часть 4. Оценка экономического эффекта модели (2 балла)
#
# В данной части мы займемся тем, что от вас скорее всего потребуется на реальной работе. А именно: мы будем считать некоторые метрики и с их помощью попытаемся настроить модель на максимизацию _прибыли_. Разумеется, здесь будет сделано множество упрощающих жизнь допущений, но обо всем по порядку. 
#
# __Задание 10 (1 балл).__ Допустим, работники вашего колл-центра получают за один звонок клиенту 1 доллар. При согласии клиента на предлагаемые условия он принесет в банк 10 долларов.
#
#
# === Краткий курс экономики от ФКН 👒===
#
#    - Если вы всё прослушали на экономике, то напомним, что выручка — это сколько денег нам принесли клиенты, а прибыль — выручка за вычетом расходов на зарплату и прочее.
#
# === Конец краткого курса экономики от ФКН 👒 ===
#
# Загрузите данные о телемаркетинге из предыдущего блока заданий. В этой части не нужно делить выборку - мы будем использовать кросс-валидацию. Используйте 5 фолдов, сделайте `shuffle=True, random_state=500`. По кросс-валидации у вас получится 5 вариантов обучающей и тестовой выборки. Обучите логистическую регрессию на каждой обучающей выборке (воспользуйтесь one-hot для категориальных признаков, гиперпараметры оставьте со значениями по умолчанию) и сделайте предсказания для соответствующих тестовых выборок. Допустим, всем положительным прогнозам ваши сотрудники решили позвонить. Посчитайте на всех тестовых выборках выручку и усредните. Сколько денег вы в среднем заработаете? Также вычислите стандартное отклонение.
#
# Сколько из заработанных денег придётся отдать операторам вашего колл-центра?

# +
#╰( ͡° ͜ʖ ͡° )つ──☆*:・ﾟ
# -

# Внесем некоторую долю случайности. Пусть теперь согласный на условия клиент будет приносить не 10 долларов, а случайную величину, равномерно распределенную в интервале $[0;20)$. Проделайте все те же самые действия. Для имитации реальной ситуации **НЕ** фиксируйте `random_seed` при подсчете выручки с клиента. Что получилось?

# +
#╰( ͡° ͜ʖ ͡° )つ──☆*:・ﾟ
# -

# Настройте по кросс-валидации коэффициент регуляризации модели для максимизации прибыли (считайте как случайную величину выше). Удалось ли получить какой-то выигрыш? При каком коэффициенте регуляризациии прибыль максимальна? Постройте график зависимости ожидаемой прибыли от коэффициента, также укажите стандартные отклонения (вам поможет `plt.errorbar`).

# +
#╰( ͡° ͜ʖ ͡° )つ──☆*:・ﾟ
# -

# __Задание 11 (1 балл).__ Выше мы уже описали примерную экономическую модель вашей задачи. Как вы считаете, что для этого бизнеса важнее — хороший precision или recall модели? Почему?
#
# __Ответ:__
#
#

# Вспомним, что на самом деле логистическая регрессия предсказывает вероятности положительного класса для объекта. Возможно, путем настройки порога бинаризации этих вероятностей мы сможем получить какой-то выигрыш? Проверьте ваши рассуждения выше с помощью настройки порога бинаризации на кросс-валидации для максимизации прибыли. Воспользуйтесь сеткой от 0 до 1 с шагом 0.01. Напомним, что снижение порога дает нам более высокий recall и более низкий precision, и наоборот.

# +
#╰( ͡° ͜ʖ ͡° )つ──☆*:・ﾟ
# -

# Каковы значения precision и recall на выбранном пороге? Оцените по кросс-валидации. Также вычислите стандартное отклонение.

# +
#╰( ͡° ͜ʖ ͡° )つ──☆*:・ﾟ
# -

# Как вы, вероятно, уже поняли, ваша модель склоняется к более высокому recall. Попробуйте оценить качество модели с помощью `PR-AUC` в зоне recall $\geq$ 0.5. Сделайте это следующим образом - выберите только те пороги, на которых достигается необходимый recall, затем интерпретируйте отсеченный в единичном квадрате прямоугольник как новый единичный квадрат и посчитайте площадь под отсеченной кривой. 

# +
#╰( ͡° ͜ʖ ͡° )つ──☆*:・ﾟ
# -

# __Бонус (1 балл):__ чтобы получить 1 балл, вставьте что угодно в ячейку ниже. Даже можно не ходить в музей.
#
# (Бонус может получить только тот, кто решил хотя бы одно задание).

#╰( ͡° ͜ʖ ͡° )つ──☆*:・ﾟ
print("!!!")

    

StarcoderdataJupyterScriptsDedupFiltered