import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline

"""
Explanation: Практическое задание к уроку 1 (2 неделя).
Линейная регрессия: переобучение и регуляризация
В этом задании мы на примерах увидим, как переобучаются линейные модели, разберем, почему так происходит, и выясним, как диагностировать и контролировать переобучение.
Во всех ячейках, где написан комментарий с инструкциями, нужно написать код, выполняющий эти инструкции. Остальные ячейки с кодом (без комментариев) нужно просто выполнить. Кроме того, в задании требуется отвечать на вопросы; ответы нужно вписывать после выделенного слова "Ответ:".
Напоминаем, что посмотреть справку любого метода или функции (узнать, какие у нее аргументы и что она делает) можно с помощью комбинации Shift+Tab. Нажатие Tab после имени объекта и точки позволяет посмотреть, какие методы и переменные есть у этого объекта.
End of explanation
"""


# (0 баллов)
# Считайте данные и выведите первые 5 строк
df = pd.read_csv('bikes_rent.csv')
df.head(5)

"""
Explanation: Мы будем работать с датасетом "bikes_rent.csv", в котором по дням записаны календарная информация и погодные условия, характеризующие автоматизированные пункты проката велосипедов, а также число прокатов в этот день. Последнее мы будем предсказывать; таким образом, мы будем решать задачу регрессии.
Знакомство с данными
Загрузите датасет с помощью функции pandas.read_csv в переменную df. Выведите первые 5 строчек, чтобы убедиться в корректном считывании данных:
End of explanation
"""


fig, axes = plt.subplots(nrows=3, ncols=4, figsize=(15, 10))
for idx, feature in enumerate(df.columns[:-1]):
    df.plot(feature, "cnt", subplots=True, kind="scatter", ax=axes[idx / 4, idx % 4])

"""
Explanation: Для каждого дня проката известны следующие признаки (как они были указаны в источнике данных):
* season: 1 - весна, 2 - лето, 3 - осень, 4 - зима
* yr: 0 - 2011, 1 - 2012
* mnth: от 1 до 12
* holiday: 0 - нет праздника, 1 - есть праздник
* weekday: от 0 до 6
* workingday: 0 - нерабочий день, 1 - рабочий день
* weathersit: оценка благоприятности погоды от 1 (чистый, ясный день) до 4 (ливень, туман)
* temp: температура в Цельсиях
* atemp: температура по ощущениям в Цельсиях
* hum: влажность
* windspeed(mph): скорость ветра в милях в час
* windspeed(ms): скорость ветра в метрах в секунду
* cnt: количество арендованных велосипедов (это целевой признак, его мы будем предсказывать)
Итак, у нас есть вещественные, бинарные и номинальные (порядковые) признаки, и со всеми из них можно работать как с вещественными. С номинальныеми признаками тоже можно работать как с вещественными, потому что на них задан порядок. Давайте посмотрим на графиках, как целевой признак зависит от остальных
End of explanation
"""


# Код 1.1 (0.5 балла)
# Посчитайте корреляции всех признаков, кроме последнего, с последним с помощью метода corrwith:
df[df.columns[:-1]].corrwith(df["cnt"])


"""
Explanation: Блок 1. Ответьте на вопросы (каждый 0.5 балла):
1. Каков характер зависимости числа прокатов от месяца? 
   * ответ: Похожа на - квадратичную, то есть к лету достигает максимума, постепенно увеличивая значения, потом после лета постепенно падает
1. Укажите один или два признака, от которых число прокатов скорее всего зависит линейно
   * ответ: от temp и от atemp
Давайте более строго оценим уровень линейной зависимости между признаками и целевой переменной. Хорошей мерой линейной зависимости между двумя векторами является корреляция Пирсона. В pandas ее можно посчитать с помощью двух методов датафрейма: corr и corrwith. Метод df.corr вычисляет матрицу корреляций всех признаков из датафрейма. Методу df.corrwith нужно подать еще один датафрейм в качестве аргумента, и тогда он посчитает попарные корреляции между признаками из df и этого датафрейма.
End of explanation
"""


# Код 1.2 (0.5 балла)
# Посчитайте попарные корреляции между признаками temp, atemp, hum, windspeed(mph), windspeed(ms) и cnt
# с помощью метода corr:
df[['temp','atemp', 'hum','windspeed(mph)','windspeed(ms)','cnt']].corr()

"""
Explanation: В выборке есть признаки, коррелирующие с целевым, а значит, задачу можно решать линейными методами.
По графикам видно, что некоторые признаки похожи друг на друга. Поэтому давайте также посчитаем корреляции между вещественными признаками.
End of explanation
"""


# Код 1.3 (0.5 балла)
# Выведите средние признаков
print df.mean()


"""
Explanation: На диагоналях, как и полагается, стоят единицы. Однако в матрице имеются еще две пары сильно коррелирующих столбцов: temp и atemp (коррелируют по своей природе) и два windspeed (потому что это просто перевод одних единиц в другие). Далее мы увидим, что этот факт негативно сказывается на обучении линейной модели.
Напоследок посмотрим средние признаков (метод mean), чтобы оценить масштаб признаков и доли 1 у бинарных признаков.
End of explanation
"""


from sklearn.preprocessing import scale
from sklearn.utils import shuffle

df_shuffled = shuffle(df, random_state=123)
X = scale(df_shuffled[df_shuffled.columns[:-1]])
y = df_shuffled["cnt"]

"""
Explanation: Признаки имеют разный масштаб, значит для дальнейшей работы нам лучше нормировать матрицу объекты-признаки.
Проблема первая: коллинеарные признаки
Итак, в наших данных один признак дублирует другой, и есть еще два очень похожих. Конечно, мы могли бы сразу удалить дубликаты, но давайте посмотрим, как бы происходило обучение модели, если бы мы не заметили эту проблему. 
Для начала проведем масштабирование, или стандартизацию признаков: из каждого признака вычтем его среднее и поделим на стандартное отклонение. Это можно сделать с помощью метода scale.
Кроме того, нужно перемешать выборку, это потребуется для кросс-валидации.
End of explanation
"""


from sklearn.linear_model import LinearRegression

# Код 2.1 (1 балл)
# Создайте объект линейного регрессора, обучите его на всех данных и выведите веса модели 
# (веса хранятся в переменной coef_ класса регрессора).
# Можно выводить пары (название признака, вес), воспользовавшись функцией zip, встроенной в язык python
# Названия признаков хранятся в переменной df.columns
alg_lnr = LinearRegression()
alg_lnr.fit(X, y)
print zip(df.columns[0:(df.shape[1]-1)], alg_lnr.coef_)

"""
Explanation: Давайте обучим линейную регрессию на наших данных и посмотрим на веса признаков.
End of explanation
"""


from sklearn.linear_model import Lasso, Ridge

# Код 2.2 (0.5 балла)
# Обучите линейную модель с L1-регуляризацией
alg_lnr = Lasso()
alg_lnr.fit(X, y)

print zip(df.columns[0:(df.shape[1]-1)], alg_lnr.coef_)

# Код 2.3 (0.5 балла)
# Обучите линейную модель с L2-регуляризацией
alg_lnr = Ridge()
alg_lnr.fit(X, y)
print zip(df.columns[0:(df.shape[1]-1)], alg_lnr.coef_)

"""
Explanation: Мы видим, что веса при линейно-зависимых признаках по модулю значительно больше, чем при других признаках.
Чтобы понять, почему так произошло, вспомним аналитическую формулу, по которой вычисляются веса линейной модели в методе наименьших квадратов:
$w = (X^TX)^{-1} X^T y$.
Если в X есть коллинеарные (линейно-зависимые) столбцы, матрица $X^TX$ становится вырожденной, и формула перестает быть корректной. Чем более зависимы признаки, тем меньше определитель этой матрицы и тем хуже аппроксимация $Xw \approx y$. Такая ситуацию называют проблемой мультиколлинеарности, вы обсуждали ее на лекции.
С парой temp-atemp чуть менее коррелирующих переменных такого не произошло, однако на практике всегда стоит внимательно следить за коэффициентами при похожих признаках.
Решение проблемы мультиколлинеарности состоит в регуляризации линейной модели. К оптимизируемому функционалу прибавляют L1 или L2 норму весов, умноженную на коэффициент регуляризации $\alpha$. В первом случае метод называется Lasso, а во втором --- Ridge. Подробнее об этом также рассказано в лекции.
Обучите регрессоры Ridge и Lasso с параметрами по умолчанию и убедитесь, что проблема с весами решилась.
End of explanation
"""


# Код 3.1 (1 балл)
alphas = np.arange(1, 500, 50)
coefs_lasso = np.zeros((alphas.shape[0], X.shape[1])) # матрица весов размера (число регрессоров) x (число признаков)
coefs_ridge = np.zeros((alphas.shape[0], X.shape[1]))
# Для каждого значения коэффициента из alphas обучите регрессор Lasso
# и запишите веса в соответствующую строку матрицы coefs_lasso (вспомните встроенную в python функцию enumerate),
# а затем обучите Ridge и запишите веса в coefs_ridge.
print coefs_lasso.shape
num = 0
for alpha  in alphas:
    modelL1 = Lasso(alpha = alpha).fit(X, y)
    modelL2 = Ridge(alpha = alpha).fit(X, y)

    coefs_lasso[num] = modelL1.coef_
    coefs_ridge[num] = modelL2.coef_
    
    num = num + 1

"""
Explanation: Блок 2. Поясните, каким образом введение регуляризации решает проблему с весами и мультиколлинеарностью.
Ваш ответ (1 балл): регуляризация изменяет функцию ошибку, таким образом модель "наказывается" за слишком большие  коэф.. Коэф. перед скоростью там получился просто коллосальный на много порядков больше, чем перед остальными, регуляризация же решила эту проблему
Проблема вторая: неинформативные признаки
В отличие от L2-регуляризации, L1 обнуляет веса при некоторых признаках. Объяснение данному факту дается в одной из лекций курса.
Давайте пронаблюдаем, как меняются веса при увеличении коэффициента регуляризации $\alpha$ (в лекции коэффициент при регуляризаторе мог быть обозначен другой буквой).
End of explanation
"""


plt.figure(figsize=(8, 5))
for coef, feature in zip(coefs_lasso.T, df.columns):
    plt.plot(alphas, coef, label=feature, color=np.random.rand(3))
plt.legend(loc="upper right", bbox_to_anchor=(1.4, 0.95))
plt.xlabel("alpha")
plt.ylabel("feature weight")
plt.title("Lasso")

plt.figure(figsize=(8, 5))
for coef, feature in zip(coefs_ridge.T, df.columns):
    plt.plot(alphas, coef, label=feature, color=np.random.rand(3))
plt.legend(loc="upper right", bbox_to_anchor=(1.4, 0.95))
plt.xlabel("alpha")
plt.ylabel("feature weight")
plt.title("Ridge")

plt.figure(figsize=(8, 5))
for coef, feature in zip(coefs_lasso.T, df.columns):
    if feature == 'windspeed(mph)' or feature == 'windspeed(ms)' or feature == 'yr':
        plt.plot(alphas, coef, label=feature, color=np.random.rand(3))
plt.legend(loc="upper right", bbox_to_anchor=(1.4, 0.95))
plt.xlabel("alpha")
plt.ylabel("feature weight")
plt.title("Lasso - temp")

"""
Explanation: Визуализируем динамику весов при увеличении параметра регуляризации:
End of explanation
"""


from sklearn.linear_model import LassoCV

# Код 3.2 (1 балл)
# Обучите регрессор LassoCV на всех параметрах регуляризации из alpha
# Постройте график _усредненного_ по строкам MSE в зависимости от alpha. 
# Выведите выбранное alpha, а также пары "признак-коэффициент" для обученного вектора коэффициентов
alphas = np.arange(1, 100, 5)

clf = LassoCV(alphas = alphas)
model = clf.fit(X, y)
print model.mse_path_.shape
plt.plot(model.alphas_, model.mse_path_.mean(axis =1))
plt.xlabel('alpha')
plt.ylabel('mean MSE')
plt.show()


print 'выбранное alpha = ', model.alpha_,'\n', zip(df.columns[0:(df.shape[1]-1)], clf.coef_)

"""
Explanation: Ответы на следующие вопросы можно давать, глядя на графики или выводя коэффициенты на печать.
Блок 3. Ответьте на вопросы (каждый 0.25 балла):
1. Какой регуляризатор (Ridge или Lasso) агрессивнее уменьшает веса при одном и том же alpha?
    * Ответ:Lasso
1. Что произойдет с весами Lasso, если alpha сделать очень большим? Поясните, почему так происходит.
    * Ответ: многие из них станут = 0, такое происходит в силу особенности функционала ошибки, при котором на ошибку влияет модуль сумму коэф.
1. Можно ли утверждать, что Lasso исключает один из признаков windspeed при любом значении alpha > 0? А Ridge? Ситается, что регуляризатор исключает признак, если коэффициент при нем < 1e-3.
    * Ответ: Да можно, это видно на 3ем графике выше.
1. Какой из регуляризаторов подойдет для отбора неинформативных признаков?
    * Ответ: L1, это можно определить по 0-ым коэф. перед ними
Далее будем работать с Lasso.
Итак, мы видим, что при изменении alpha модель по-разному подбирает коэффициенты признаков. Нам нужно выбрать наилучшее alpha. 
Для этого, во-первых, нам нужна метрика качества. Будем использовать в качестве метрики сам оптимизируемый функционал метода наименьших квадратов, то есть Mean Square Error.
Во-вторых, нужно понять, на каких данных эту метрику считать. Нельзя выбирать alpha по значению MSE на обучающей выборке, потому что тогда мы не сможем оценить, как модель будет делать предсказания на новых для нее данных. Если мы выберем одно разбиение выборки на обучающую и тестовую (это называется holdout), то настроимся на конкретные "новые" данные, и вновь можем переобучиться. Поэтому будем делать несколько разбиений выборки, на каждом пробовать разные значения alpha, а затем усреднять MSE. Удобнее всего делать такие разбиения кросс-валидацией, то есть разделить выборку на K частей, или блоков, и каждый раз брать одну из них как тестовую, а из оставшихся блоков составлять обучающую выборку. 
Делать кросс-валидацию для регрессии в sklearn совсем просто: для этого есть специальный регрессор, LassoCV, который берет на вход список из alpha и для каждого из них вычисляет MSE на кросс-валидации. После обучения (если оставить параметр cv=3 по умолчанию) регрессор будет содержать переменную mse_path_, матрицу размера len(alpha) x k, k = 3 (число блоков в кросс-валидации), содержащую значения MSE на тесте для соответствующих запусков. Кроме того, в переменной alpha_ будет храниться выбранное значение параметра регуляризации, а в coef_, традиционно, обученные веса, соответствующие этому alpha_.
Обратите внимание, что регрессор может менять порядок, в котором он проходит по alphas; для сопоставления с матрицей MSE лучше использовать переменную регрессора alphas_.
End of explanation
"""


# Код 3.3 (1 балл)
# Выведите значения alpha, соответствующие минимумам MSE на каждом разбиении (то есть по столбцам).
# На трех отдельных графиках визуализируйте столбцы .mse_path_

print model.alphas_[np.argmin(model.mse_path_[:,0])]
plt.figure(figsize=(8, 5))
plt.plot(model.alphas_, model.mse_path_[:,0], 'r')
plt.xlabel('alpha')
plt.ylabel('1st MSE')
plt.legend()
plt.axhline(y= np.amin(model.mse_path_[:,0],axis =0))
plt.title('Min alpha = ' + str(model.alphas_[np.argmin(model.mse_path_[:,0])]))

print model.alphas_[np.argmin(model.mse_path_[:,1])]

plt.figure(figsize=(8, 5))
plt.plot(model.alphas_, model.mse_path_[:,1], color = 'y')
plt.xlabel('alpha')
plt.ylabel('2nd MSE')
plt.axhline(y= np.amin(model.mse_path_[:,1],axis =0))
plt.title('Min alpha = ' + str(model.alphas_[np.argmin(model.mse_path_[:,1])]))

print model.alphas_[np.argmin(model.mse_path_[:,2])]

plt.figure(figsize=(8, 5))
plt.plot(model.alphas_, model.mse_path_[:,2], color = 'g')
plt.xlabel('alpha')
plt.ylabel('3rd MSE')
plt.axhline(y= np.amin(model.mse_path_[:,2],axis =0))
plt.title('Min alpha = ' + str(model.alphas_[np.argmin(model.mse_path_[:,2])]))



"""
Explanation: Итак, мы выбрали некоторый параметр регуляризации. Давайте посмотрим, какие бы мы выбирали alpha, если бы делили выборку только один раз на обучающую и тестовую, то есть рассмотрим траектории MSE, соответствующие отдельным блокам выборки.
End of explanation
"""


zip(df.columns[0:(df.shape[1]-1)], clf.coef_)

"""
Explanation: На каждом разбиении оптимальное значение alpha свое, и ему соответствует большое MSE на других разбиениях. Получается, что мы настраиваемся на конкретные обучающие и контрольные выборки. При выборе alpha на кросс-валидации мы выбираем нечто "среднее", что будет давать приемлемое значение метрики на разных разбиениях выборки. 
Наконец, как принято в анализе данных, давайте проинтерпретируем результат.
Блок 4. Ответьте на вопросы (каждый 0.5 балла):
1. В последней обученной модели выберите 4 признака с наибольшими (положительными) коэфициентами (и выпишите их), посмотрите на визуализации зависимостей cnt от этих признаков, которые мы рисовали в блоке "Знакомство с данными". Видна ли возрастающая линейная зависимость cnt от этих признаков по графикам? Логично ли утверждать (из здравого смысла), что чем больше значение этих признаков, тем больше людей захотят взять велосипеды? 
    * Ответ: yr, season, atemp, temp. Чем больше значение признака, тем больше cnt, но почему так - не так очевидно. Рост температуры и проката связан с тем, что просто зимой не покатаешься, летом же это уже не так важно. Осенью меньше чем летом прокат. А второй год больше значение в силу роста самого бизнеса. Не все так однозначно.
1. Выберите 3 признака с наибольшими по модулю отрицательными коэффициентами (и выпишите их), посмотрите на соответствующие визуализации. Видна ли убывающая линейная зависимость? Логично ли утверждать, что чем больше величина этих признаков, тем меньше людей захотят взять велосипеды?
    * Ответ: weathersit, hum, windspeed. Да видна - очевидно при дожде и ветре кататься не хочется, как и при плохом прогнозе на погоду.
1. Выпишите признаки с коэффициентами, близкими к нулю (< 1e-3). Как вы думаете, почему модель исключила их из модели (вновь посмотрите на графики)? Верно ли, что они никак не влияют на спрос на велосипеды?
    * Ответ: windspeed(ms) - этот признак влияет, просто его лассо обрезало регуляризацией в силу мультиколлинеарности, что уже обсуждалось.
End of explanation
"""


