<jupyter_start><jupyter_text>Введение в анализ данных, ИАД-1 НИУ ВШЭ, 2018-19 учебный год  Домашнее задание №2 Задание выполнил(а): Подчезерцев <NAME>  Общая информация __Дата выдачи:__ 25.02.2019__Дедлайн:__ 23:59 5.03.2019  О задании В данном домашнем задании вы реализуете линейную регрессию своими руками и сравните её с версией в scikit-learn.  Оценивание и штрафыЗа сдачу задания позже срока на итоговую оценку за задание накладывается штраф в размере 1 балл в день, но получить отрицательную оценку нельзя. __Внимание!__ Домашнее задание выполняется самостоятельно. «Похожие» решения считаются плагиатом и все задействованные студенты (в том числе те, у кого списали) не могут получить за него больше 0 баллов.  Формат сдачи __Стирать условия нельзя!__Загрузка файлов с решениями происходит в системе [Anytask](https://anytask.org/).Формат названия файла: homework_02_Фамилия_Имя.ipynb  Часть 1. Обыкновенная линейная регрессия<jupyter_code>import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

%matplotlib inline<jupyter_output><empty_output><jupyter_text>1 (4 балла). На семинаре мы начали реализовывать класс линейной регрессии; теперь необходимо его закончить.<jupyter_code>class LinReg():
    def __init__(self, step_l=0.017, step_p=1.3, step_s=29, num_steps=100, eps=1e-6):
        self.w = None # веса
        self.w_curr = None
        self.shape = None
        self.num_steps = num_steps
        self.eps = eps
        self.step_l = step_l
        self.step_p = step_p
        self.step_s = step_s
    
    def _calc_grad(self, X_train, y_train):
        return 2 * np.dot(X_train.T, np.dot(X_train, self.w) - y_train) / self.shape[1]
    
    def _calc_step(self, step):
        return self.step_l*(self.step_s/(self.step_s+step))**self.step_p
    
    def _calc_weight(self, X_train, y_train, step):
        return self.w - self._calc_step(step) * self._calc_grad(X_train, y_train)
    
    def fit(self, X_train, y_train):
        # self.w = np.random.random(self.count)
        self.shape = X_train.shape
        self.w = np.zeros(self.shape[1])
        for step in range(self.num_steps):
            self.w_curr = self._calc_weight(X_train, y_train, step)
            if np.linalg.norm(self.w_curr - self.w) < self.eps:
                break
            self.w = self.w_curr
            
        return self
        
    def predict(self, X_test):
        return np.dot(X_test, self.w)<jupyter_output><empty_output><jupyter_text>2 (3 балла). Проверим корректность работы класса на датасете Boston Housing. * Загрузите его из sklearn* Проведите предобработку данных(нормализация вещественных признаков, one-hot encoding категориальных)* Разделите на обучение и контроль в соотношении 80:20 с random_seed 42 (самый популярный в прошлом задании)* Обучите собственную реализацию регрессии и сравните по метрикам MAE и RMSE с LinearRegression из sklearn.<jupyter_code>from sklearn.datasets import load_boston
dataset = load_boston()
X = pd.DataFrame(dataset['data'], columns=dataset['feature_names'])
y = dataset['target']
X["MEDV"] = y
X["ONES"] = 1.0
def norm(data, col):
    data[col] = (data[col] - data[col].mean())/data[col].std()
X.head()
norm(X, "CRIM")
norm(X, "ZN")
norm(X, "INDUS")
norm(X, "NOX")
norm(X, "RM")
norm(X, "AGE")
norm(X, "DIS")
norm(X, "RAD")
norm(X, "TAX")
norm(X, "PTRATIO")
norm(X, "B")
norm(X, "LSTAT")
norm(X, "MEDV")
y = X["MEDV"]
del X["MEDV"]
X.head()
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, mean_absolute_error
def RMSE(test, predict):
    return np.sqrt(mean_squared_error(test, predict))
lr = LinearRegression()
lr.fit(X_train,y_train)
predict = lr.predict(X_test)
print('RMSE: %f' % RMSE(y_test, predict))
print('MAE: %f' % mean_absolute_error(y_test, predict))
sk_w = lr.coef_
sk_w
lr = LinReg(0.014,1,29)
lr.fit(X_train,y_train)
predict = lr.predict(X_test)
print('RMSE: %f' % RMSE(y_test, predict))
print('MAE: %f' % mean_absolute_error(y_test, predict))
my_w = lr.w
my_w
print("Compate weights")
print('RMSE: %f' % RMSE(sk_w, my_w))
print('MAE: %f' % mean_absolute_error(sk_w, my_w))<jupyter_output>Compate weights
RMSE: 0.011758
MAE: 0.007314<jupyter_text>Величины RMSE и MAE практически совпали для выборок. Метрики для весов отличаются на сотые доли  3 (3 балла). Линейная регрессия зачастую легко переобучается - модель необходимо штрафовать за величину весов; для этого применяют L1 и L2 регуляризацию: добавление нормы весов к функции потерь. В случае L2-регулязации функционал будет выглядеть как$$ L = (Xw - y)^T(Xw - y) + \lambda||w||_2 $$.Параметр $\lambda$ подбирается на отложенной выборке или по кросс-валидации.* Реализуйте обучение линейной регрессии с L2-регуляризацией* Найдите оптимальный с точки зрения метрики MAE коэффициент $\lambda$ (по кросс-валидации)* Постройте график зависимости метрики на тестовой выборке от $\lambda$ (__подпишите оси__)* Сравните результаты с Ridge регрессией из sklearn (аналогично пункту 2).<jupyter_code>class RegL2(LinReg):
    def __init__(self, lamb, step_l=0.017, step_p=1.3, step_s=29, num_steps=100, eps=1e-6):
        LinReg.__init__(self, step_l, step_p, step_s, num_steps, eps)
        self.lamb = lamb

    def _calc_grad(self, X_train, y_train):
        return 2 * np.dot(X_train.T, np.dot(X_train, self.w) - y_train) / self.shape[1] + 2 * self.lamb * self.w
from sklearn.linear_model import Ridge
from sklearn.model_selection import KFold
def iterator(start, end, step):
    i = start
    while i <= end:
        yield i
        i += step
for n in range(2,11):
    plt_dataX = []
    plt_dataY = []
    kf = KFold(n_splits=n,shuffle=True, random_state=42)
    for i in iterator(0, 6, 0.1):
        plt_dataX.append(i)
        err = 0
        for train_index, test_index in kf.split(X):
            X_train, X_test = X.values[train_index], X.values[test_index]
            y_train, y_test = y[train_index], y[test_index]
            lr_l2 = RegL2(i)
            lr_l2.fit(X_train,y_train)
            err += mean_absolute_error(y_test, lr_l2.predict(X_test))
        plt_dataY.append(err / kf.get_n_splits())
    plt.figure(figsize=(20,10))
    plt.grid(True)
    plt.title('Качество модели от коэффициента гамма, N=' + str(n))
    plt.xlabel('Гамма')
    plt.ylabel('MAE')
    plt.scatter(plt_dataX, plt_dataY)
    plt.show()
    print("Best lambda value:", plt_dataX[plt_dataY.index(min(plt_dataY))], "MAE:", min(plt_dataY))<jupyter_output><empty_output><jupyter_text>Начиная с некоторого значения разбиения N стабильность зависимости MAE от гаммы падаетРазобьем выборку на $\sqrt{size(X)}$ частей и возьмем данные в качестве основной<jupyter_code>kf = KFold(n_splits=round(X.shape[1]**0.5),shuffle=True, random_state=42)
plt_dataX = []
plt_dataY = []
for i in iterator(0, 10, 0.1):
    plt_dataX.append(i)
    err = 0
    for train_index, test_index in kf.split(X):
        X_train, X_test = X.values[train_index], X.values[test_index]
        y_train, y_test = y[train_index], y[test_index]
        lr_l2 = RegL2(i)
        lr_l2.fit(X_train,y_train)
        err += mean_absolute_error(y_test, lr_l2.predict(X_test))
    plt_dataY.append(err / kf.get_n_splits())
plt.figure(figsize=(20,10))
plt.grid(True)
plt.title('Качество модели от коэффициента гамма, N=' + str(kf.get_n_splits()))
plt.xlabel('Гамма')
plt.ylabel('MAE')
plt.scatter(plt_dataX, plt_dataY)
plt.show()
print("Best lambda value:", plt_dataX[plt_dataY.index(min(plt_dataY))], "MAE:", min(plt_dataY))<jupyter_output><empty_output><jupyter_text>Результаты на стандартном разбиении 2 задания<jupyter_code>X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
plt_dataX = []
plt_dataY = []
for i in iterator(0, 10, 0.1):
    plt_dataX.append(i)
    lr_l2 = RegL2(i)
    lr_l2.fit(X_train,y_train)
    plt_dataY.append(mean_absolute_error(y_test, lr_l2.predict(X_test)))
plt.figure(figsize=(20,10))
plt.grid(True)
plt.title('Качество модели от коэффициента гамма')
plt.xlabel('Гамма')
plt.ylabel('MAE')
plt.scatter(plt_dataX, plt_dataY)
plt.show()
print("Best lambda value:", plt_dataX[plt_dataY.index(min(plt_dataY))], "MAE:", min(plt_dataY))
lr_l2 = RegL2(2.3)
lr_l2.fit(X_train,y_train)
predict = lr_l2.predict(X_test)
print('RMSE: %f' % RMSE(y_test, predict))
print('MAE: %f' % mean_absolute_error(y_test, predict))
my_w = lr_l2.w
my_w
lr_l2 = Ridge()
lr_l2.fit(X_train,y_train)
predict = lr_l2.predict(X_test)
print('RMSE: %f' % RMSE(y_test, predict))
print('MAE: %f' % mean_absolute_error(y_test, predict))
sk_w = lr_l2.coef_
sk_w
print("Compate weights")
print('RMSE: %f' % RMSE(sk_w, my_w))
print('MAE: %f' % mean_absolute_error(sk_w, my_w))<jupyter_output>Compate weights
RMSE: 0.065430
MAE: 0.047945<jupyter_text>Ошибки почти совпали, веса отличаются на сотые доли  4 (Бонус, 4 балла). * Реализуйте обучение линейной регрессии с L1-регуляризацией (класс RegL1)* Проведите сравнения аналогично пункту 3 (но с Lasso вместо Ridge).<jupyter_code>from sklearn.linear_model import Lasso
class RegL1(LinReg):
    def __init__(self, lamb, step_l=0.017, step_p=1.3, step_s=29, num_steps=100, eps=1e-6):
        LinReg.__init__(self, step_l, step_p, step_s, num_steps, eps)
        self.lamb = lamb

    def _calc_grad(self, X_train, y_train):
        return 2 * np.dot(X_train.T, np.dot(X_train, self.w) - y_train) / self.shape[1] + self.lamb * np.sign(self.w)
kf = KFold(n_splits=round(X.shape[1]**0.5),shuffle=True, random_state=42)
plt_dataX = []
plt_dataY = []
for i in iterator(0, 2, 0.01):
    plt_dataX.append(i)
    err = 0
    for train_index, test_index in kf.split(X):
        X_train, X_test = X.values[train_index], X.values[test_index]
        y_train, y_test = y[train_index], y[test_index]
        lr_l1 = RegL1(i)
        lr_l1.fit(X_train,y_train)
        err += mean_absolute_error(y_test, lr_l1.predict(X_test))
    plt_dataY.append(err / kf.get_n_splits())
plt.figure(figsize=(20,10))
plt.grid(True)
plt.title('Качество модели от коэффициента лямбда, N=' + str(kf.get_n_splits()))
plt.xlabel('лямбда')
plt.ylabel('MAE')
plt.scatter(plt_dataX, plt_dataY)
plt.show()
print("Best lambda value:", plt_dataX[plt_dataY.index(min(plt_dataY))], "MAE:", min(plt_dataY))
kf = KFold(n_splits=round(X.shape[1]**0.5),shuffle=True, random_state=42)
plt_dataX = []
plt_dataY = []
for i in iterator(0, 0.02, 0.001):
    plt_dataX.append(i)
    err = 0
    for train_index, test_index in kf.split(X):
        X_train, X_test = X.values[train_index], X.values[test_index]
        y_train, y_test = y[train_index], y[test_index]
        lr_l1 = Lasso(i)
        lr_l1.fit(X_train,y_train)
        err += mean_absolute_error(y_test, lr_l1.predict(X_test))
    plt_dataY.append(err / kf.get_n_splits())
plt.figure(figsize=(20,10))
plt.grid(True)
plt.title('Качество модели от коэффициента лямбда, N=' + str(kf.get_n_splits()))
plt.xlabel('лямбда')
plt.ylabel('MAE')
plt.scatter(plt_dataX, plt_dataY)
plt.show()
print("Best lambda value:", plt_dataX[plt_dataY.index(min(plt_dataY))], "MAE:", min(plt_dataY))
lr_l1 = Lasso(0.008)
lr_l1.fit(X_train,y_train)
predict = lr_l1.predict(X_test)
print('RMSE: %f' % RMSE(y_test, predict))
print('MAE: %f' % mean_absolute_error(y_test, predict))
sk_w = lr_l1.coef_
sk_w
lr_l1 = RegL1(0.8)
lr_l1.fit(X_train,y_train)
predict = lr_l1.predict(X_test)
print('RMSE: %f' % RMSE(y_test, predict))
print('MAE: %f' % mean_absolute_error(y_test, predict))
my_w = lr_l1.w
my_w
print("Compate weights")
print('RMSE: %f' % RMSE(sk_w, my_w))
print('MAE: %f' % mean_absolute_error(sk_w, my_w))<jupyter_output>Compate weights
RMSE: 0.012384
MAE: 0.008634<jupyter_text>Модели почти похожи, коэффициенты и ошибки совпадают  5 (Бонус, 2 балла). Исследуйте для реализации регрессии с L2-регуляризацией зависимость качества на тестовой выборке (с графиками) от:* Длины шага* Количества шагов спуска* Константы epsilon.<jupyter_code>class RegL2Spec(LinReg):
    def __init__(self, lamb=2.3, step=0.001, num_steps=100, eps=1e-6):
        LinReg.__init__(self, None, None, None, num_steps, eps)
        self.lamb = lamb
        self.step = step
        
    def _calc_step(self, step):
        return self.step
    
    def _calc_grad(self, X_train, y_train):
        return 2 * np.dot(X_train.T, np.dot(X_train, self.w) - y_train) / self.shape[1] + 2 * self.lamb * self.w
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
plt_dataX = []
plt_dataY = []
for i in iterator(10, 100, 1):
    plt_dataX.append(i)
    lr_l2 = RegL2Spec(num_steps=i)
    lr_l2.fit(X_train,y_train)
    plt_dataY.append(mean_absolute_error(y_test, lr_l2.predict(X_test)))
for i in iterator(100, 1000, 10):
    plt_dataX.append(i)
    lr_l2 = RegL2Spec(num_steps=i)
    lr_l2.fit(X_train,y_train)
    plt_dataY.append(mean_absolute_error(y_test, lr_l2.predict(X_test)))
plt.figure(figsize=(20,10))
plt.grid(True)
plt.title('Качество модели от количества шагов')
plt.xlabel('Количества шагов')
plt.ylabel('MAE')
plt.scatter(plt_dataX, plt_dataY)
plt.show()
print("Best value:", plt_dataX[plt_dataY.index(min(plt_dataY))], "MAE:", min(plt_dataY))<jupyter_output><empty_output><jupyter_text>Улучшение модели снижается около 100 итерации, практически перестает изменятся после 300-400<jupyter_code>plt_dataX = []
plt_dataY = []
for i in iterator(0.0001, 0.006, 0.0001):
    lr_l2 = RegL2Spec(num_steps=100, step = i)
    lr_l2.fit(X_train,y_train)
    if mean_absolute_error(y_test, lr_l2.predict(X_test)) < 1:
        plt_dataX.append(i)
        plt_dataY.append(mean_absolute_error(y_test, lr_l2.predict(X_test)))
plt.figure(figsize=(20,10))
plt.grid(True)
plt.title('Качество модели от размера шага')
plt.xlabel('Шаг')
plt.ylabel('MAE')
plt.scatter(plt_dataX, plt_dataY)
plt.show()
print("Best value:", plt_dataX[plt_dataY.index(min(plt_dataY))], "MAE:", min(plt_dataY))<jupyter_output><empty_output><jupyter_text>Валидный размер шага практически не влияет на качество модели.На границах интервала наблюдается резкое ухудшенеи качество модели. Для малого шага обосновывается недостаточным временем для обучения, для больших - выход из области минимума.Наиболее оптимальный один из самых большых шагов<jupyter_code>plt_dataX = []
plt_dataY = []
def log_iterator(start, end, step):
    i = start
    while i <= end:
        yield i
        i *= step
for i in log_iterator(10**-7, 1, 1.4):
    lr_l2 = RegL2Spec(num_steps=100, step = 0.0054, eps=i)
    lr_l2.fit(X_train,y_train)
    plt_dataX.append(i)
    plt_dataY.append(mean_absolute_error(y_test, lr_l2.predict(X_test)))
plt.figure(figsize=(20,10))
plt.grid(True)
plt.title('Качество модели от эпсилон')
plt.xlabel('Eps')
plt.ylabel('MAE')
plt.scatter(plt_dataX, plt_dataY)
plt.xscale("log", nonposx='clip')
plt.xlim(10**-7)
plt.show()
print("Best value:", plt_dataX[plt_dataY.index(min(plt_dataY))], "MAE:", min(plt_dataY))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1324383891003807, 'fraction_numerical': 0.02604688439190543, 'mean_word_length': 4.717449408171058, 'pattern_counts': {'":': 0, '<': 43, '<?xml version=': 0, '>': 39, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '313215', 'license': 'Apache-2.0', 'n_tokens_mistral': 6679, 'n_tokens_neox': 6752, 'n_words': 1253, 'path': 'introduction_to_machine_learning/hw_02.ipynb', 'repo_name': 'AsciiShell/hse_ida', 'size': '14985'}