Задание основано на материалах лекций по линейной регрессии и градиентному спуску. Вы будете прогнозировать выручку компании в зависимости от уровня ее инвестиций в рекламу по TV, в газетах и по радио.
Линейная регрессия - один из наиболее хорошо изученных методов машинного обучения, позволяющий прогнозировать значения количественного признака в виде линейной комбинации прочих признаков с параметрами - весами модели. Оптимальные (в смысле минимальности некоторого функционала ошибки) параметры линейной регрессии можно найти аналитически с помощью нормального уравнения или численно с помощью методов оптимизации.
Линейная регрессия использует простой функционал качества - среднеквадратичную ошибку. Мы будем работать с выборкой, содержащей 3 признака. Для настройки параметров (весов) модели решается следующая задача: $$\Large \frac{1}{\ell}\sum_{i=1}^\ell{{((w_0 + w_1x_{i1} + w_2x_{i2} + w_3x_{i3}) - y_i)}^2} \rightarrow \min_{w_0, w_1, w_2, w_3},$$ где $x_{i1}, x_{i2}, x_{i3}$ - значения признаков $i$-го объекта, $y_i$ - значение целевого признака $i$-го объекта, $\ell$ - число объектов в обучающей выборке.
Параметры $w_0, w_1, w_2, w_3$, по которым минимизируется среднеквадратичная ошибка, можно находить численно с помощью градиентного спуска. Градиентный шаг для весов будет выглядеть следующим образом: $$\Large w_0 \leftarrow w_0 - \frac{2\eta}{\ell} \sum_{i=1}^\ell{{((w_0 + w_1x_{i1} + w_2x_{i2} + w_3x_{i3}) - y_i)}}$$ $$\Large w_j \leftarrow w_j - \frac{2\eta}{\ell} \sum_{i=1}^\ell{{x_{ij}((w_0 + w_1x_{i1} + w_2x_{i2} + w_3x_{i3}) - y_i)}},\ j \in \{1,2,3\}$$ Здесь $\eta$ - параметр, шаг градиентного спуска.
Проблема градиентного спуска, описанного выше, в том, что на больших выборках считать на каждом шаге градиент по всем имеющимся данным может быть очень вычислительно сложно. В стохастическом варианте градиентного спуска поправки для весов вычисляются только с учетом одного случайно взятого объекта обучающей выборки: $$\Large w_0 \leftarrow w_0 - \frac{2\eta}{\ell} {((w_0 + w_1x_{k1} + w_2x_{k2} + w_3x_{k3}) - y_k)}$$ $$\Large w_j \leftarrow w_j - \frac{2\eta}{\ell} {x_{kj}((w_0 + w_1x_{k1} + w_2x_{k2} + w_3x_{k3}) - y_k)},\ j \in \{1,2,3\},$$ где $k$ - случайный индекс, $k \in \{1, \ldots, \ell\}$.
Нахождение вектора оптимальных весов $w$ может быть сделано и аналитически. Мы хотим найти такой вектор весов $w$, чтобы вектор $y$, приближающий целевой признак, получался умножением матрицы $X$ (состоящей из всех признаков объектов обучающей выборки, кроме целевого) на вектор весов $w$. То есть, чтобы выполнялось матричное уравнение: $$\Large y = Xw$$ Домножением слева на $X^T$ получаем: $$\Large X^Ty = X^TXw$$ Это хорошо, поскольку теперь матрица $X^TX$ - квадратная, и можно найти решение (вектор $w$) в виде: $$\Large w = {(X^TX)}^{-1}X^Ty$$ Матрица ${(X^TX)}^{-1}X^T$ - псевдообратная для матрицы $X$. В NumPy такую матрицу можно вычислить с помощью функции numpy.linalg.pinv.
Однако, нахождение псевдообратной матрицы - операция вычислительно сложная и нестабильная в случае малого определителя матрицы $X$ (проблема мультиколлинеарности). На практике лучше находить вектор весов $w$ решением матричного уравнения $$\Large X^TXw = X^Ty$$Это может быть сделано с помощью функции numpy.linalg.solve.
Но все же на практике для больших матриц $X$ быстрее работает градиентный спуск, особенно его стохастическая версия.
В начале напишем простую функцию для записи ответов в текстовый файл. Ответами будут числа, полученные в ходе решения этого задания, округленные до 3 знаков после запятой. Полученные файлы после выполнения задания надо отправить в форму на странице задания на Coursera.org.
In [2]:
def write_answer_to_file(answer, filename):
with open(filename, 'w') as f_out:
f_out.write(str(round(answer, 3)))
1. Загрузите данные из файла advertising.csv в объект pandas DataFrame. Источник данных.
In [3]:
import pandas as pd
adver_data = pd.read_csv('advertising.csv')
Посмотрите на первые 5 записей и на статистику признаков в этом наборе данных.
In [4]:
adver_data.head(5)
Out[4]:
In [5]:
adver_data.describe()
Out[5]:
Создайте массивы NumPy X из столбцов TV, Radio и Newspaper и y - из столбца Sales. Используйте атрибут values объекта pandas DataFrame.
In [6]:
import numpy as np
X = adver_data[['TV', 'Radio', 'Newspaper']].values
y = adver_data[['Sales']].values
print 'X array', X.shape, ':'
print X
print 'y array', y.shape, ':'
print y
Отмасштабируйте столбцы матрицы X, вычтя из каждого значения среднее по соответствующему столбцу и поделив результат на стандартное отклонение. Для определенности, используйте методы mean и std векторов NumPy (реализация std в Pandas может отличаться). Обратите внимание, что в numpy вызов функции .mean() без параметров возвращает среднее по всем элементам массива, а не по столбцам, как в pandas. Чтобы произвести вычисление по столбцам, необходимо указать параметр axis.
In [7]:
means, stds = X.mean(axis=0), X.std(axis=0)
print 'Mean values:', means
print 'Std. deviations:', stds
In [8]:
X = (X - means) / stds
print 'X array', X.shape, ':'
print X
Добавьте к матрице X столбец из единиц, используя методы hstack, ones и reshape библиотеки NumPy. Вектор из единиц нужен для того, чтобы не обрабатывать отдельно коэффициент $w_0$ линейной регрессии.
In [9]:
X = np.hstack( (X, np.ones((X.shape[0], 1))) )
print 'X array', X.shape, ':'
print X
2. Реализуйте функцию mserror - среднеквадратичную ошибку прогноза. Она принимает два аргумента - объекты Series y (значения целевого признака) и y_pred (предсказанные значения). Не используйте в этой функции циклы - тогда она будет вычислительно неэффективной.
In [10]:
def mserror(y, y_pred):
return np.sum((y - y_pred) ** 2) / y.shape[0]
Какова среднеквадратичная ошибка прогноза значений Sales, если всегда предсказывать медианное значение Sales по исходной выборке? Запишите ответ в файл '1.txt'.
In [11]:
answer1 = mserror(y, np.median(y))
print 'Median sales value:', np.median(y)
print(answer1)
write_answer_to_file(answer1, '1.txt')
3. Реализуйте функцию normal_equation, которая по заданным матрицам (массивам NumPy) X и y вычисляет вектор весов $w$ согласно нормальному уравнению линейной регрессии.
In [12]:
def normal_equation(X, y):
return np.dot(np.linalg.pinv(X), y)
In [13]:
norm_eq_weights = normal_equation(X, y)
print(norm_eq_weights)
Какие продажи предсказываются линейной моделью с весами, найденными с помощью нормального уравнения, в случае средних инвестиций в рекламу по ТВ, радио и в газетах? (то есть при нулевых значениях масштабированных признаков TV, Radio и Newspaper). Запишите ответ в файл '2.txt'.
In [14]:
answer2 = norm_eq_weights[3]
print(answer2)
write_answer_to_file(answer2, '2.txt')
4. Напишите функцию linear_prediction, которая принимает на вход матрицу X и вектор весов линейной модели w, а возвращает вектор прогнозов в виде линейной комбинации столбцов матрицы X с весами w.
In [15]:
def linear_prediction(X, w):
return np.dot(X, w)
Какова среднеквадратичная ошибка прогноза значений Sales в виде линейной модели с весами, найденными с помощью нормального уравнения? Запишите ответ в файл '3.txt'.
In [16]:
answer3 = mserror(y, linear_prediction(X, norm_eq_weights))
print(answer3)
write_answer_to_file(answer3, '3.txt')
5. Напишите функцию stochastic_gradient_step, реализующую шаг стохастического градиентного спуска для линейной регрессии. Функция должна принимать матрицу X, вектора y и w, число train_ind - индекс объекта обучающей выборки (строки матрицы X), по которому считается изменение весов, а также число $\eta$ (eta) - шаг градиентного спуска (по умолчанию eta=0.01). Результатом будет вектор обновленных весов. Наша реализация функции будет явно написана для данных с 3 признаками, но несложно модифицировать для любого числа признаков, можете это сделать.
In [76]:
def stochastic_gradient_step(X, y, w, train_ind, eta=0.01):
#массив градиента для каждого параметра
grad = np.zeros( (X.shape[1], 1) )
#перевод строки параметров в столбец для поэлементного умножения
xRowToCol = X[train_ind].reshape(-1, 1)
#вычисление градиента для каждого параметра
for i in xrange(X.shape[1]):
grad[i] = xRowToCol[i] * ( np.sum(xRowToCol * w) - y[train_ind])
return (2 * eta / X.shape[0]) * grad
6. Напишите функцию stochastic_gradient_descent, реализующую стохастический градиентный спуск для линейной регрессии. Функция принимает на вход следующие аргументы:
На каждой итерации в вектор (список) должно записываться текущее значение среднеквадратичной ошибки. Функция должна возвращать вектор весов $w$, а также вектор (список) ошибок.
In [77]:
def stochastic_gradient_descent(X, y, w_init, eta=1e-2, max_iter=1e4,
max_weight_dist=1e-8, seed=42, verbose=False):
# Инициализируем расстояние между векторами весов на соседних
# итерациях большим числом.
weight_dist = np.inf
# Инициализируем вектор весов
w = w_init
# Сюда будем записывать ошибки на каждой итерации
errors = []
# Счетчик итераций
iter_num = 0
# Будем порождать псевдослучайные числа
# (номер объекта, который будет менять веса), а для воспроизводимости
# этой последовательности псевдослучайных чисел используем seed.
np.random.seed(seed)
# Основной цикл
while weight_dist > max_weight_dist and iter_num < max_iter:
# порождаем псевдослучайный
# индекс объекта обучающей выборки
random_ind = np.random.randint(X.shape[0])
#вычисляем значение весов на следующей итерации
wOld = w
w = wOld - stochastic_gradient_step(X, y, w, random_ind, eta=eta)
#расстояние между векторами весов на соседних итерациях
weight_dist = np.linalg.norm(wOld - w)
#среднеквадратичная ошибка для найденного вектора весов
err = mserror( y, linear_prediction(X, w) )
errors.append(err)
#счётчик итераций
iter_num += 1
return w, errors
Запустите $10^5$ итераций стохастического градиентного спуска. Укажите вектор начальных весов w_init, состоящий из нулей. Оставьте параметры eta и seed равными их значениям по умолчанию (eta=0.01, seed=42 - это важно для проверки ответов).
In [78]:
%%time
stoch_grad_desc_weights, stoch_errors_by_iter = stochastic_gradient_descent(X, y, np.zeros( (X.shape[1], 1) ), max_iter=1e5)
Посмотрим, чему равна ошибка на первых 50 итерациях стохастического градиентного спуска. Видим, что ошибка не обязательно уменьшается на каждой итерации.
In [88]:
%pylab inline
plot(range(200), stoch_errors_by_iter[:200])
xlabel('Iteration number')
ylabel('MSE')
Out[88]:
Теперь посмотрим на зависимость ошибки от номера итерации для $10^5$ итераций стохастического градиентного спуска. Видим, что алгоритм сходится.
In [80]:
%pylab inline
plot(range(len(stoch_errors_by_iter)), stoch_errors_by_iter)
axis([0, 1e5, 0, 250])
grid(True)
xlabel('Iteration number')
ylabel('MSE')
Out[80]:
Посмотрим на вектор весов, к которому сошелся метод.
In [81]:
stoch_grad_desc_weights
Out[81]:
Посмотрим на среднеквадратичную ошибку на последней итерации.
In [82]:
stoch_errors_by_iter[-1]
Out[82]:
Какова среднеквадратичная ошибка прогноза значений Sales в виде линейной модели с весами, найденными с помощью градиентного спуска? Запишите ответ в файл '4.txt'.
In [83]:
answer4 = mserror(y, linear_prediction(X, stoch_grad_desc_weights))
print(answer4)
write_answer_to_file(answer4, '4.txt')
Ответами к заданию будут текстовые файлы, полученные в ходе этого решения. Обратите внимание, что отправленные файлы не должны содержать пустую строку в конце. Данный нюанс является ограничением платформы Coursera. Мы работаем над исправлением этого ограничения.