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