Линейная регрессия и основные библиотеки Python для анализа данных и научных вычислений

Это задание посвящено линейной регрессии. На примере прогнозирования роста человека по его весу Вы увидите, какая математика за этим стоит, а заодно познакомитесь с основными библиотеками Python, необходимыми для дальнейшего прохождения курса.

Материалы

Задание 1. Первичный анализ данных c Pandas

В этом заданиии мы будем использовать данные SOCR по росту и весу 25 тысяч подростков.

[1]. Если у Вас не установлена библиотека Seaborn - выполните в терминале команду conda install seaborn. (Seaborn не входит в сборку Anaconda, но эта библиотека предоставляет удобную высокоуровневую функциональность для визуализации данных).


In [1]:
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import scipy.optimize as opt
%matplotlib inline

Считаем данные по росту и весу (weights_heights.csv, приложенный в задании) в объект Pandas DataFrame:


In [2]:
data = pd.read_csv('weights_heights.csv', index_col='Index')

Чаще всего первое, что надо надо сделать после считывания данных - это посмотреть на первые несколько записей. Так можно отловить ошибки чтения данных (например, если вместо 10 столбцов получился один, в названии которого 9 точек с запятой). Также это позволяет познакомиться с данными, как минимум, посмотреть на признаки и их природу (количественный, категориальный и т.д.).

После этого стоит построить гистограммы распределения признаков - это опять-таки позволяет понять природу признака (степенное у него распределение, или нормальное, или какое-то еще). Также благодаря гистограмме можно найти какие-то значения, сильно не похожие на другие - "выбросы" в данных. Гистограммы удобно строить методом plot Pandas DataFrame с аргументом kind='hist'.

Пример. Построим гистограмму распределения роста подростков из выборки data. Используем метод plot для DataFrame data c аргументами y='Height' (это тот признак, распределение которого мы строим)


In [3]:
data.plot(y='Height', kind='hist', 
           color='red',  title='Height (inch.) distribution')


Out[3]:
<matplotlib.axes._subplots.AxesSubplot at 0xc5ea278>

Аргументы:

  • y='Height' - тот признак, распределение которого мы строим
  • kind='hist' - означает, что строится гистограмма
  • color='red' - цвет

[2]. Посмотрите на первые 5 записей с помощью метода head Pandas DataFrame. Нарисуйте гистограмму распределения веса с помощью метода plot Pandas DataFrame. Сделайте гистограмму зеленой, подпишите картинку.


In [4]:
# Ваш код здесь
# Первые 5 записей из данных
data.head(5)


Out[4]:
Height Weight
Index
1 65.78331 112.9925
2 71.51521 136.4873
3 69.39874 153.0269
4 68.21660 142.3354
5 67.78781 144.2971

In [5]:
# Ваш код здесь
# Диаграмма расределения веса
data.plot(y='Weight', kind='hist', 
           color='green',  title='Weight (inch.) distribution')


Out[5]:
<matplotlib.axes._subplots.AxesSubplot at 0xca3e438>

Один из эффективных методов первичного анализа данных - отображение попарных зависимостей признаков. Создается $m \times m$ графиков (m - число признаков), где по диагонали рисуются гистограммы распределения признаков, а вне диагонали - scatter plots зависимости двух признаков. Это можно делать с помощью метода $scatter\_matrix$ Pandas Data Frame или pairplot библиотеки Seaborn.

Чтобы проиллюстрировать этот метод, интересней добавить третий признак. Создадим признак Индекс массы тела (BMI). Для этого воспользуемся удобной связкой метода apply Pandas DataFrame и lambda-функций Python.


In [6]:
def make_bmi(height_inch, weight_pound):
    METER_TO_INCH, KILO_TO_POUND = 39.37, 2.20462
    return (weight_pound / KILO_TO_POUND) / \
           (height_inch / METER_TO_INCH) ** 2

In [7]:
data['BMI'] = data.apply(lambda row: make_bmi(row['Height'], 
                                              row['Weight']), axis=1)

[3]. Постройте картинку, на которой будут отображены попарные зависимости признаков , 'Height', 'Weight' и 'BMI' друг от друга. Используйте метод pairplot библиотеки Seaborn.


In [8]:
# Ваш код здесь
g = sns.pairplot(data)#,diag_kind="hist",kind="reg"


Часто при первичном анализе данных надо исследовать зависимость какого-то количественного признака от категориального (скажем, зарплаты от пола сотрудника). В этом помогут "ящики с усами" - boxplots библиотеки Seaborn. Box plot - это компактный способ показать статистики вещественного признака (среднее и квартили) по разным значениям категориального признака. Также помогает отслеживать "выбросы" - наблюдения, в которых значение данного вещественного признака сильно отличается от других.

[4]. Создайте в DataFrame data новый признак weight_category, который будет иметь 3 значения: 1 – если вес меньше 120 фунтов. (~ 54 кг.), 3 - если вес больше или равен 150 фунтов (~68 кг.), 2 – в остальных случаях. Постройте «ящик с усами» (boxplot), демонстрирующий зависимость роста от весовой категории. Используйте метод boxplot библиотеки Seaborn и метод apply Pandas DataFrame. Подпишите ось y меткой «Рост», ось x – меткой «Весовая категория».


In [9]:
# определение значения для признака weight_category
def weight_category(weight):
    if weight<120:
        return 1
    elif weight>=150:
        return 3
    else:
        return 2

data['weight_cat'] = data['Weight'].apply(weight_category)
# Ваш код здесь
sns.boxplot(x=data["weight_cat"], y=data["Height"])
data.head(5)


Out[9]:
Height Weight BMI weight_cat
Index
1 65.78331 112.9925 18.357573 1
2 71.51521 136.4873 18.762577 2
3 69.39874 153.0269 22.338895 3
4 68.21660 142.3354 21.504526 2
5 67.78781 144.2971 22.077581 2

[5]. Постройте scatter plot зависимости роста от веса, используя метод plot для Pandas DataFrame с аргументом kind='scatter'. Подпишите картинку.


In [10]:
# Ваш код здесь
data.plot(x='Weight',y='Height', kind='scatter')
plt.title('Correlation Weight and Height')


Out[10]:
<matplotlib.text.Text at 0xeb84710>

Задание 2. Минимизация квадратичной ошибки

В простейшей постановке задача прогноза значения вещественного признака по прочим признакам (задача восстановления регрессии) решается минимизацией квадратичной функции ошибки.

[6]. Напишите функцию, которая по двум параметрам $w_0$ и $w_1$ вычисляет квадратичную ошибку приближения зависимости роста $y$ от веса $x$ прямой линией $y = w_0 + w_1 * x$: $$error(w_0, w_1) = \sum_{i=1}^n {(y_i - (w_0 + w_1 * x_i))}^2 $$ Здесь $n$ – число наблюдений в наборе данных, $y_i$ и $x_i$ – рост и вес $i$-ого человека в наборе данных.


In [32]:
# Ваш код здесь
# ошибка от 2-х параметров
def error(w1,w0,data):
    sum = 0.0
    for i in xrange(1,len(data)+1):
        sum += (data.loc[i,"Height"] - (w0 + w1 * data.loc[i,"Weight"]))**2
    return sum

Итак, мы решаем задачу: как через облако точек, соответсвующих наблюдениям в нашем наборе данных, в пространстве признаков "Рост" и "Вес" провести прямую линию так, чтобы минимизировать функционал из п. 6. Для начала давайте отобразим хоть какие-то прямые и убедимся, что они плохо передают зависимость роста от веса.

[7]. Проведите на графике из п. 5 Задания 1 две прямые, соответствующие значениям параметров ($w_0, w_1) = (60, 0.05)$ и ($w_0, w_1) = (50, 0.16)$. Используйте метод plot из matplotlib.pyplot, а также метод linspace библиотеки NumPy. Подпишите оси и график.


In [34]:
# Ваш код здесь
x = np.linspace(70,180)

def y(w0,w1,x):
    return w0 + w1 * x

data.plot(x='Weight',y='Height', kind='scatter')
plt.title('Correlation Weight and Height')
plt.plot(x,y(60,0.05,x),color='r')
plt.plot(x,y(50,0.16,x),color='g')
plt.legend( ('w0 = 60, w1 = 0.05', 'w0 = 50, w1 = 0.16') )


Out[34]:
<matplotlib.legend.Legend at 0x1067cc18>

Минимизация квадратичной функции ошибки - относительная простая задача, поскольку функция выпуклая. Для такой задачи существует много методов оптимизации. Посмотрим, как функция ошибки зависит от одного параметра (наклон прямой), если второй параметр (свободный член) зафиксировать.

[8]. Постройте график зависимости функции ошибки, посчитанной в п. 6, от параметра $w_1$ при $w_0$ = 50. Подпишите оси и график.


In [35]:
# Ваш код здесь
# ошибка от параметра w1
w0 = 50
err_graph = []
for w1 in range(-10,10,1):
    err_graph.append(error(w1,w0,data))
    plt.plot(err_graph)


Теперь методом оптимизации найдем "оптимальный" наклон прямой, приближающей зависимость роста от веса, при фиксированном коэффициенте $w_0 = 50$.

[9]. С помощью метода minimize_scalar из scipy.optimize найдите минимум функции, определенной в п. 6, для значений параметра $w_1$ в диапазоне [-5,5]. Проведите на графике из п. 5 Задания 1 прямую, соответствующую значениям параметров ($w_0$, $w_1$) = (50, $w_1\_opt$), где $w_1\_opt$ – найденное в п. 8 оптимальное значение параметра $w_1$.


In [36]:
# Ваш код здесь
w0 = 50
w1_opt = opt.minimize_scalar(error,bounds=(-5, 5),args=(w0,data))
w1_opt.x


Out[36]:
0.14109203728834385

In [43]:
# Ваш код здесь
x = np.linspace(70,180)

def y(w0,w1,x):
    return w0 + w1 * x

data.plot(y='Height',x='Weight', kind='scatter')
plt.plot(x,y(w0,w1_opt.x,x),color='r')


Out[43]:
[<matplotlib.lines.Line2D at 0x110a5048>]

При анализе многомерных данных человек часто хочет получить интуитивное представление о природе данных с помощью визуализации. Увы, при числе признаков больше 3 такие картинки нарисовать невозможно. На практике для визуализации данных в 2D и 3D в данных выделаяют 2 или, соответственно, 3 главные компоненты (как именно это делается - мы увидим далее в курсе) и отображают данные на плоскости или в объеме.

Посмотрим, как в Python рисовать 3D картинки, на примере отображения функции $z(x,y) = sin(\sqrt{x^2+y^2})$ для значений $x$ и $y$ из интервала [-5,5] c шагом 0.25.


In [45]:
from mpl_toolkits.mplot3d import Axes3D

Создаем объекты типа matplotlib.figure.Figure (рисунок) и matplotlib.axes._subplots.Axes3DSubplot (ось).


In [46]:
fig = plt.figure()
ax = fig.gca(projection='3d') # get current axis

# Создаем массивы NumPy с координатами точек по осям X и У. 
# Используем метод meshgrid, при котором по векторам координат 
# создается матрица координат. Задаем нужную функцию Z(x, y).
X = np.arange(-5, 5, 0.25)
Y = np.arange(-5, 5, 0.25)
X, Y = np.meshgrid(X, Y)
Z = np.sin(np.sqrt(X**2 + Y**2))

# Наконец, используем метод *plot_surface* объекта 
# типа Axes3DSubplot. Также подписываем оси.
surf = ax.plot_surface(X, Y, Z)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()


[10]. Постройте 3D-график зависимости функции ошибки, посчитанной в п.6 от параметров $w_0$ и $w_1$. Подпишите ось $x$ меткой «Intercept», ось $y$ – меткой «Slope», a ось $z$ – меткой «Error».


In [49]:
fig = plt.figure()
ax = fig.gca(projection='3d') # get current axis

# Создаем массивы NumPy с координатами точек по осям X и У. 
# Используем метод meshgrid, при котором по векторам координат 
# создается матрица координат. Задаем нужную функцию Z(x, y).
X = np.arange(-10,10,1)#w0
Y = np.arange(-10,10,1)#w1
X, Y = np.meshgrid(X, Y)
Z = error(X,Y,data)

# Наконец, используем метод *plot_surface* объекта 
# типа Axes3DSubplot. Также подписываем оси.
surf = ax.plot_surface(X, Y, Z)
ax.set_xlabel('Intercept')
ax.set_ylabel('Slope')
ax.set_zlabel('Error')
ax.set_title('Error in 3D')
plt.show()



In [19]:
# Ваш код здесь

In [20]:
# Ваш код здесь

[11]. С помощью метода minimize из scipy.optimize найдите минимум функции, определенной в п. 6, для значений параметра $w_0$ в диапазоне [-100,100] и $w_1$ - в диапазоне [-5, 5]. Начальная точка – ($w_0$, $w_1$) = (0, 0). Используйте метод оптимизации L-BFGS-B (аргумент method метода minimize). Проведите на графике из п. 5 Задания 1 прямую, соответствующую найденным оптимальным значениям параметров $w_0$ и $w_1$. Подпишите оси и график.


In [58]:
# Ваш код здесь
def func_err(w_opt,data):
    return error(w_opt[0],w_opt[1],data)

w_opt = opt.minimize(func_err,[0.0,0.0],method = 'L-BFGS-B',
                     bounds=([-5,5],[-100,100]),args=(data))
w_opt


Out[58]:
      fun: 67545.287085287142
 hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>
      jac: array([ 0.18335413,  0.00727596])
  message: 'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
     nfev: 51
      nit: 12
   status: 0
  success: True
        x: array([  0.08200696,  57.57171564])

In [61]:
# Ваш код здесь
x = np.linspace(70,180)

def y(w0,w1,x):
    return w0 + w1 * x

data.plot(y='Height',x='Weight', kind='scatter')
plt.plot(x,y(w_opt.x[1],w_opt.x[0],x),color='r')


Out[61]:
[<matplotlib.lines.Line2D at 0x106ad588>]

Критерии оценки работы

  • Выполняется ли тетрадка IPython без ошибок? (15 баллов)
  • Верно ли отображена гистограмма распределения роста из п. 2? (3 балла). Правильно ли оформлены подписи? (1 балл)
  • Верно ли отображены попарные зависимости признаков из п. 3? (3 балла). Правильно ли оформлены подписи? (1 балл)
  • Верно ли отображена зависимость роста от весовой категории из п. 4? (3 балла). Правильно ли оформлены подписи? (1 балл)
  • Верно ли отображен scatter plot роста от веса из п. 5? (3 балла). Правильно ли оформлены подписи? (1 балл)
  • Правильно ли реализована функция подсчета квадратичной ошибки из п. 6? (10 баллов)
  • Правильно ли нарисован график из п. 7? (3 балла) Правильно ли оформлены подписи? (1 балл)
  • Правильно ли нарисован график из п. 8? (3 балла) Правильно ли оформлены подписи? (1 балл)
  • Правильно ли используется метод minimize_scalar из scipy.optimize? (6 баллов). Правильно ли нарисован график из п. 9? (3 балла) Правильно ли оформлены подписи? (1 балл)
  • Правильно ли нарисован 3D-график из п. 10? (6 баллов) Правильно ли оформлены подписи? (1 балл)
  • Правильно ли используется метод minimize из scipy.optimize? (6 баллов). Правильно ли нарисован график из п. 11? (3 балла). Правильно ли оформлены подписи? (1 балл)