#!/usr/bin/env python
# coding: utf-8
# # Линейная регрессия
# Шестаков А.В. Майнор по анализу данных - 01/03/2016
# Сегодня мы рассмотрим следующие темы:
#
# 1. Постановка задачи линейной регрессии
# 2. Преобразование переменных и интерпретация
# 3. Переобучение\недообучение, мультиколлинеарность и регуляризация
# ## 1. Постановка задачи
# Дано описание $n$ объектов по $m$ признакам. Обычно оно выражается в виде матрицы размера $n \times m$: $X = [x^{(i)}_j]^{i=1\dots n}_{j=1\dots m} $.
($x^{(i)}_j$ означает $j$-ый признак $i$-го объекта)
# Дана **вещественная** зависимая переменная, которая тоже имеет отношение к этим объекам: $y \in \mathbb{R}^n$.
#
# Наша задача, выявить **линейную** зависимость между признаками в $X$ и значениями в $y$:
# $$\hat{y} = X\beta \quad \Leftrightarrow \quad \hat{y}^{(i)} = \beta_0 + \beta_1x^{(i)}_1 + \dots$$
# То есть необходимо оценить коэффициенты $\beta_i$.
#
# В случае линейной регрессии коэффициенты $\beta_i$ рассчитываются так, чтобы минимизировать сумму квадратов ошибок по всем наблюдениям:
# $$ L(\beta) = \frac{1}{2n}(\hat{y} - y)^{\top}(\hat{y} - y) = \frac{1}{2n}(X\beta - y)^{\top}(X\beta - y) \rightarrow \min$$ $$ \Updownarrow $$ $$ L(\beta_0,\beta_1,\dots) = \frac{1}{2n}\sum^{n}_{i=1}(\hat{y}^{(i)} - y^{(i)})^2 = \frac{1}{2n}\sum^{n}_{i=1}(\beta_0 + \beta_1x^{(i)}_1 + \dots - y^{(i)})^2 \rightarrow \min $$
#
# Мы уже знаем несколько способов решения этой задачи (Семинар 4):
# * Градиентный спуск
# * Normal Equations (Проекционные матрицы)
# ### Библиотеки для расчетов
# Самыми распространенными библиотеками в `Python` для работы с линейными методами (регрессии и классификации) являются [`scikit-learn`](http://scikit-learn.org/stable/) и [`statmodels`](http://statsmodels.sourceforge.net/). Так как в дальнейшем мы скорее всего будем работать со `scikit-learn`, то и примеры по большей части будут демонстрироваться именно в нем.
# ### Пример: Стоимость автомобиля
# Загрузите [данные](http://bit.ly/1gIQs6C) по характеристикам автомобилей Honda Accord. Названия столбцов говорят сами за себя.
# In[ ]:
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
plt.style.use('ggplot')
get_ipython().run_line_magic('matplotlib', 'inline')
# In[ ]:
df = pd.read_csv('http://bit.ly/1gIQs6C')
# In[ ]:
df.head()
# Выберем одну переменную mileage в качестве предиктора, а переменную price в качестве зависимой переменной
# In[ ]:
y = df.price.values
X = df.mileage.values.reshape(-1,1)
# В последних версиях sklearn начинает ругаться на одномерные данные (когда array.shape = (m,))
# Поэтому с помошщью .reshape(-1,1) мы искусственно добавляем единичную размерность к Х
# In[ ]:
from sklearn.linear_model import LinearRegression
# In[ ]:
model = LinearRegression(fit_intercept=True)
model.fit(X, y)
# In[ ]:
print 'Модель:\nprice = %.2f + (%.2f)*mileage' % (model.intercept_, model.coef_[0])
# In[ ]:
x = np.linspace(0, max(X), 100)
y_line = model.intercept_ + model.coef_[0]*x
fig, ax = plt.subplots(1, 1, figsize=(10,5))
ax.scatter(X, y)
ax.plot(x, y_line, c='red')
# ### Остатки и меры качества
# Давайте взглянем на ошибки (остатки)
# In[ ]:
fig, ax = plt.subplots(1, 2, figsize=(14, 5))
y_hat = model.predict(X)
res = y - y_hat
ax[0].hist(res)
ax[0].set_xlabel('residuals')
ax[0].set_ylabel('counts')
ax[1].scatter(X, res)
ax[1].set_xlabel('mileage')
ax[1].set_ylabel('residuals')
# Важно смотреть на остатки.
# Во-первых, они должны быть нормально распределены (Теорема Гаусса-Маркова).
# Во-вторых, не должно быть ярких зависимостей между значениями признака и остатками.
#
# Посмотрим на меры качества
# In[ ]:
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import median_absolute_error
from sklearn.metrics import mean_squared_error
from sklearn.metrics import r2_score
# Можно посчитать простые варианты агрегирования остатков, например:
#
# * $\frac{1}{n} \sum_i |\hat{y}^{(i)}-y^{(i)}|$ - средняя абсолютная ошибка
# * $\frac{1}{n} \sum_i (\hat{y}^{(i)}-y^{(i)})^2$ - средняя квадратичная ошибка
# In[ ]:
print 'Средняя абсолютная ошибка %.2f' % mean_absolute_error(y, y_hat)
print 'Средняя квадратичная ошибка %.2f' % mean_squared_error(y, y_hat)
# Можно рассмотреть более сложную меру: коэффициент детерминации $R^2$:
#
# * $TSS = \sum_i (y^{(i)}-\bar{y})^2$ - общая сумма квадратов (total sum of squares)
# * $RSS = \sum_i (\hat{y}^{(i)}-y^{(i)})^2$ - сумма квадратов остатков (residual sum of squares)
# * $ESS = \sum_i (\hat{y}^{(i)}-\bar{y})^2$ - объясненная сумма квадратов (explained sum of squares)
#
# Для простоты будем считать, что
# $$TSS = ESS + RSS$$
#
# Тогда Коэффициент детерминации $R^2=1-\frac{RSS}{TSS}$
#
# In[ ]:
print 'R^2 %.2f:' % r2_score(y, y_hat)
# Тоже самое через `statmodels`
# In[ ]:
import statsmodels.api as sm
model = sm.OLS(y, sm.add_constant(X))
results = model.fit()
print(results.summary())
# ## 2. Преобразование переменных
# ### Нормализация
# Как мы с вами убедились из ДЗ2, может решить всё)
# Переход к близким или единым шкалам улучшает сходимость градиентного спуска, уменьшает риск переполнения разрядности чисел, однако приходится жертвовать прямой интерпретируемостью..
#
# Нормализацию обычно проделывают для вещественных признаков.
#
# Нормализация z-score:
# 1. Вычитаем среднее: $x - \bar{x}$
# 2. Делим на стандартное отклонение: $\frac{x - \bar{x}}{std(x)}$
#
# Можно проделать вручную, можно с помошью метода ниже
# In[ ]:
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
X = scaler.fit_transform(X)
model = LinearRegression(fit_intercept=True)
model.fit(X, y)
# In[ ]:
print 'Модель:\nprice = %.2f + (%.2f)*mileage`' % (model.intercept_, model.coef_[0])
# ### Номинальная шкала
# In[ ]:
df.head()
# Обратим внимание, что в нашем DataFrame есть и другие переменные - год, особенности внешности, количество цилиндров, тип коробки передач. Как учесть их в модели? Иными словами, как их перекодировать, для использование в модели регресии?
#
# Часть из них можно представить в виде так называемых dummy (фиктивных) переменных
# In[ ]:
from sklearn.feature_extraction import DictVectorizer
cols = df.columns[1:]
dv = DictVectorizer()
dv.fit(df[cols].T.to_dict().values())
X = dv.transform(df[cols].T.to_dict().values())
# Давайте поймем, как устроен X
# In[ ]:
# #### Пример: Стоимость автомобиля - больше переменных!
# Обратите внимание, что в случае, когда мы добавляем все категории в матрицу `X`, рассчитывать вес для свободного члена не нужно.
#
# Причиной тому - так называемая Dummy Variable Trap. Матрица $X^\top X$ перестаёт быть обратимой из-за линейной зависимости столбцов (своего рода мультиколлинеарность).
# Хотя, даже если поставить `fit_intercept=True`, sklearn вас даже не остановит...
# In[ ]:
model = LinearRegression(fit_intercept=False)
model.fit(X, y)
# In[ ]:
# Выведите значение коэффициента для каждой переменной и проинтерпретируйте их значение
# In[ ]:
y_hat = model.predict(X)
print 'Средняя абсолютная ошибка %.2f' % mean_absolute_error(y, y_hat)
print 'Средняя квадратичная ошибка %.2f' % mean_squared_error(y, y_hat)
print 'R^2 %.2f:' % r2_score(y, y_hat)
# ### Природа зависимости
# Далеко не всегда переменные зависят друг от друга именно в том виде, в котором они даны. Никто не запрещает зависимость вида
# $$\log(y) = \beta_0 + \beta_1\log(x_1)$$
# или
# $$y = \beta_0 + \beta_1\frac{1}{x_1}$$
# или
# $$y = \beta_0 + \beta_1\log(x_1)$$
# или
# $$y = \beta_0 + \beta_1 x_1^2 + \beta_2 x_2^2 + \beta_3 x_1x_2 $$
# и т.д.
#
# Не смотря на то, что могут возникать какие-то нелинейные функции - всё это сводится к **линейной** регрессии (например, о втором пункте, произведите замену $z_1 = \frac{1}{x_1}$)
# #### Пример: Вес тела - мозгов
# Загрузите [данные](https://www.dropbox.com/s/8srfeh34lnj2cb3/weights.csv?dl=0) и информацией о весах мозга и тел различных биологических видов. Вес тела задан в килограммах, вес могза в граммах.
# In[ ]:
df = pd.read_csv('weights.csv', sep=';', index_col=0)
df.head()
# In[ ]:
df.plot(x = 'body_w', y='brain_w', kind='scatter')
for k, v in df.iterrows():
plt.annotate(k, v[:2])
# Должно получится что-то несуразное..
# Теперь давайте возьмем логарифм от обеих переменных и сонова нарисуем их на графике
# In[ ]:
# Your Code Here
# Постройде линейную регрессию над логарифмами значений. Найдите коэффициенты и проинтерпретируйте их
# In[ ]:
# Your Code Here
# ## 3. Переобучение\недообучение, мультиколлинеарность и регуляризация
# Одна из важнейших характеристик моделей, будь то линейная регрессия, наивные Байес и др. - их **обобщающая способность**.
# Наша задача не построить "идеальную" модель, на имеющихся у нас наблюдениях, которая идеально их будет предсказывать, но и применять эту модель для новых данных.
#
# Ниже приводятся примеры 3х моделей.
#
# [Andrew's Ng Machine Learning Class - Stanford]
# Второй момент, который важен для линейных моделей - **мультиколлинеарность**. Этот эффект возникает, когда пара предикторов близка к взаимной линейной зависимости (коэффициент корреляции по модулю близок к 1). Из-за этого:
#
# * Матрица $X^{\top} X$ становится плохо обусловленной или необратимой
# * Зависимость $y = \beta_0 + \beta_1x_1 + \beta_2x_2$ перестаёт быть одназначной
#
# С этим эффектом можно бороться несколькими способами
#
# * Последовательно добавлять переменные в модель
# * Исключать коррелируемые предикторы
# ### Регуляризация
# В обоих случаях может помочь **регуляризация** - добавление штрафного слагаемого за сложность модели в функцию потерь. В случае линейной регрессии было:
# $$ L(\beta_0,\beta_1,\dots) = \frac{1}{2n}\sum^{n}_{i=1}(\hat{y}^{(i)} - y^{(i)})^2 $$
# Стало (Ridge Regularization)
# $$ L(\beta_0,\beta_1,\dots) = \frac{1}{2n}[ \sum^{n}_{i=1}(\hat{y}^{(i)} - y^{(i)})^2 + \lambda\sum_{j=1}^{m}\beta_j^2]$$
# или (Lasso Regularization)
# $$ L(\beta_0,\beta_1,\dots) = \frac{1}{2n}[ \sum^{n}_{i=1}(\hat{y}^{(i)} - y^{(i)})^2 + \lambda\sum_{j=1}^{m}|\beta_j|]$$
#
#
#
# In[ ]:
# В sklearn эти методы называются так
from sklearn.linear_model import Lasso
from sklearn.linear_model import Ridge