Базовые принципы машинного обучения на примере линейной регрессии20.02.2017 16:18

Мне выпала честь сделать первый пост, и я, пожалуй, отклонюсь от своей привычной нейросетевой тематики и сделаю пост о базовых понятиях машинного обучения на примере одной из самых простых и самых полезных моделей — линейной регрессии. Я буду использовать язык питон для демонстрации экспериментов и отрисовки графиков, все это вы с легкостью сможете повторить на своем компьютере. Поехали.

Формализмы

Машинное обучение — это подраздел искусственного интеллекта, в котором изучаются алгоритмы, способные обучаться без прямого программирования того, что нужно изучать. Линейная регрессия является типичным представителем алгоритмов машинного обучения. Для начала ответим на вопрос «а что вообще значит обучаться?». Ответ на этот вопрос мы возьмем из книги 1997 года (стоит отметить, что оглавление этой книги не сильно отличается от современных книг по машинному обучению).

Говорят, что программа обучается на опыте относительно класса задач в смысле меры качества , если при решении задачи качество, измеряемое мерой , возрастает при демонстрации нового опыта .

Можно выделить следующие задачи , решаемые машинным обучением: обучение с учителем, обучение без учителя, обучение с подкреплением, активное обучение, трансфер знаний и т.д. Регрессия (как и классификация) относится к классу задач обучения с учителем, когда по заданному набору признаков наблюдаемого объекта необходимо спрогнозировать некоторую целевую переменную. Как правило, в задачах обучения с учителем, опыт представляется в виде множества пар признаков и целевых переменных: . В случае линейной регрессии признаковое описание объекта — это действительный вектор , а целевая переменная — это скаляр . Самой простой мерой качества для задачи регрессии является , где  — это наша оценка реального значения целевой переменной.

У нас есть задача, данные и способ оценки программы/модели. Давайте определим, что такое модель, и что значит обучить модель. Предиктивная модель — это параметрическое семейство функций (семейство гипотез):

• 
•  — множество параметров

1. обучение: ;
2. применение: .

Линейная регрессия

def generate_wave_set(n_support=1000, n_train=25, std=0.3):
    data = {}
    # выберем некоторое количество точек из промежутка от 0 до 2*pi
    data['support'] = np.linspace(0, 2*np.pi, num=n_support)
    # для каждой посчитаем значение sin(x) + 1
    # это будет ground truth
    data['values'] = np.sin(data['support']) + 1
    # из support посемплируем некоторое количество точек с возвратом, это будут признаки
    data['x_train'] = np.sort(np.random.choice(data['support'], size=n_train, replace=True))
    # опять посчитаем sin(x) + 1 и добавим шум, получим целевую переменную
    data['y_train'] = np.sin(data['x_train']) + 1 + np.random.normal(0, std, size=data['x_train'].shape[0])
    return data

data = generate_wave_set(1000, 250)

print 'Shape of X is', data['x_train'].shape
print 'Head of X is', data['x_train'][:10]

margin = 0.3
plt.plot(data['support'], data['values'], 'b--', alpha=0.5, label='manifold')
plt.scatter(data['x_train'], data['y_train'], 40, 'g', 'o', alpha=0.8, label='data')
plt.xlim(data['x_train'].min() - margin, data['x_train'].max() + margin)
plt.ylim(data['y_train'].min() - margin, data['y_train'].max() + margin)
plt.legend(loc='upper right', prop={'size': 20})
plt.title('True manifold and noised data')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

А теперь реализуем алгоритм обучения, используя магию NumPy:

# добавим колонку единиц к единственному столбцу признаков
X = np.array([np.ones(data['x_train'].shape[0]), data['x_train']]).T
# перепишем, полученную выше формулу, используя numpy
# шаг обучения - в этом шаге мы ищем лучшую гипотезу h
w = np.dot(np.dot(np.linalg.inv(np.dot(X.T, X)), X.T), data['y_train'])
# шаг применения: посчитаем прогноз
y_hat = np.dot(w, X.T)

margin = 0.3
plt.plot(data['support'], data['values'], 'b--', alpha=0.5, label='manifold')
plt.scatter(data['x_train'], data['y_train'], 40, 'g', 'o', alpha=0.8, label='data')

plt.plot(data['x_train'], y_hat, 'r', alpha=0.8, label='fitted')

plt.xlim(data['x_train'].min() - margin, data['x_train'].max() + margin)
plt.ylim(data['y_train'].min() - margin, data['y_train'].max() + margin)
plt.legend(loc='upper right', prop={'size': 20})
plt.title('Fitted linear regression')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

Как мы видим, линия не очень-то совпадает с настоящей кривой. Среднеквадратичная ошибка равна 0.26704 условных единиц. Очевидно, что если бы вместо линии мы использовали кривую третьего порядка, то результат был бы куда лучше. И, на самом деле, с помощью линейной регрессии мы можем обучать нелинейные модели.

Полиномиальная регрессия

# список степеней p полиномов, который мы протестируем
degree_list = [1, 2, 3, 5, 7, 10, 13]

cmap = plt.get_cmap('jet')
colors = [cmap(i) for i in np.linspace(0, 1, len(degree_list))]

margin = 0.3
plt.plot(data['support'], data['values'], 'b--', alpha=0.5, label='manifold')
plt.scatter(data['x_train'], data['y_train'], 40, 'g', 'o', alpha=0.8, label='data')

w_list = []
err = []
for ix, degree in enumerate(degree_list):
    # список с предрасчитанными степенями признака
    dlist = [np.ones(data['x_train'].shape[0])] + \
                map(lambda n: data['x_train']**n, range(1, degree + 1))
    X = np.array(dlist).T
    w = np.dot(np.dot(np.linalg.inv(np.dot(X.T, X)), X.T), data['y_train'])
    w_list.append((degree, w))
    y_hat = np.dot(w, X.T)
    err.append(np.mean((data['y_train'] - y_hat)**2))
    plt.plot(data['x_train'], y_hat, color=colors[ix], label='poly degree: %i' % degree)

plt.xlim(data['x_train'].min() - margin, data['x_train'].max() + margin)
plt.ylim(data['y_train'].min() - margin, data['y_train'].max() + margin)
plt.legend(loc='upper right', prop={'size': 20})
plt.title('Fitted polynomial regressions')
plt.xlabel('x')
plt.ylabel('y')
plt.show() 

На графике мы можем наблюдать сразу два феномена. Пока не обращайте внимание на 13-ую степень полинома. При увеличении степени полинома, средняя ошибка продолжает уменьшаться, хотя мы вроде были уверены, что именно кубический полном должен лучше всего описывать наши данные.

Это явный признак переобучения, который можно заметить по визуализации даже не используя тестовый набор данных: при увеличении степени полинома выше третьей модель начинает интерполировать данные, вместо экстраполяции. Другими словами, график функции проходит точно через точки из тренировочного набора данных, причем чем выше степень полинома, тем через большее количество точек он проходит. Степень полинома отражает сложность модели. Таким образом, сложные модели, у которых степеней свободы достаточно много, могут попросту запомнить весь тренировочный набор, полностью теряя обобщающую способность. Это и есть проявление негативной стороны принципа минимизации эмпирического риска.

Вернемся к полиному 13-ой степени, с ним явно что-то не так. По идее, мы ожидаем, что полином 13-ой степени будет описывать тренировочный набор данных еще лучше, но результат показывает, что это не так. Из курса линейной алгебры мы помним, что обратная матрица существует только для несингулярных матриц, т.е. тех, у которых нет линейной зависимости колонок или строк. В методе наименьших квадратов нам необходимо инвертировать матрицу ковариации признаков: . Для тестирования на линейную зависимость или мультиколлинеарность можно использовать число обусловленности матрицы. Один из способов оценки этого числа для матриц — это отношение модуля максимального собственного числа матрицы к модулю минимального собственного числа. Большое число обусловленности матрицы, или же наличие одного или нескольких собственных чисел близких к нулю свидетельствует о наличии мультиколлинеарности (или нечеткой мультиколлиниарности, когда ). Такие матрицы называются слабо обусловленными, а задача — некорректно поставленной. При инвертировании такой матрицы, решения имеют большую дисперсию. Это проявляется в том, что при небольшом изменении начальной матрицы, инвертированные будут сильно отличаться друг от друга. На практике это всплывет тогда, когда к 1000 семплов, вы добавите всего один, а решение МНК будет совсем другим. Посмотрим на собственные числа матрицы ковариации, нас там ждет сюрприз:

np.linalg.eigvals(np.cov(X[:, 1:].T))

Все так, numpy вернул два комплекснозначных собственных значения, что идет вразрез с теорией. Для симметричных и положительно определенных матриц (каковой и является матрица ковариации) все собственные значения должны быть действительные. Возможно, это произошло из-за того, что при работе с большими числами матрица стала слегка несимметричной, но это не точно ¯\_(ツ)_/¯. Если вы вдруг найдете причину такого поведения нумпая, пожалуйста, напишите в комменте.

Перед тем как перейти к следующему разделу, давайте посмотрим на амплитуду параметров полиномиальной регрессии. Мы увидим, что при увеличении степени полинома, размах значений коэффициентов растет чуть ли не экспоненциально. Да, они еще и скачут в разные стороны.

for ix, t in enumerate(w_list):
    degree, w = t
    fig, ax = plt.subplots()
    plt.bar(range(max(degree_list) + 1), np.hstack((w, [0]*(max(degree_list) - w.shape[0] + 1))), color=colors[ix])
    plt.title('Magnitude of fitted LR parameters for poly:%i' % degree)
    plt.xlabel('degree')
    plt.ylabel('value of w')
    ax.set_xticks(np.array(range(max(degree_list) + 1)) + 0.5)
    ax.set_xticklabels(range(max(degree_list) + 1))
    plt.show()

Регуляризация

•  — это коэффициент регуляризации, то, насколько сильно мы хотим учитывать условие

•  — единичная диагональна матрица

Такая регрессия называется гребневой регрессией (ridge regression). А гребнем является как раз диагональная матрица которую мы прибавляем к матрице ковариации с линейнозависимыми колонками, в результате получаемая матрица не сингулярна.

Для такой матрицы число обусловленности будет равно: , где  — это собственные числа матрицы. Таким образом, увеличивая параметр регуляризации мы уменьшаем число обусловленности, а обусловленность задачи улучшается.

# define regularization parameter
lmbd = 0.1

degree_list = [1, 2, 3, 10, 12, 13]
cmap = plt.get_cmap('jet')
colors = [cmap(i) for i in np.linspace(0, 1, len(degree_list))]

margin = 0.3
plt.plot(data['support'], data['values'], 'b--', alpha=0.5, label='manifold')
plt.scatter(data['x_train'], data['y_train'], 40, 'g', 'o', alpha=0.8, label='data')

w_list_l2 = []
err = []
for ix, degree in enumerate(degree_list):
    dlist = [[1]*data['x_train'].shape[0]] + map(lambda n: data['x_train']**n, range(1, degree + 1))
    X = np.array(dlist).T
    w = np.dot(np.dot(np.linalg.inv(np.dot(X.T, X) + lmbd*np.eye(X.shape[1])), X.T), data['y_train'])
    w_list_l2.append((degree, w))
    y_hat = np.dot(w, X.T)
    plt.plot(data['x_train'], y_hat, color=colors[ix], label='poly degree: %i' % degree)
    err.append(np.mean((data['y_train'] - y_hat)**2))

plt.xlim(data['x_train'].min() - margin, data['x_train'].max() + margin)
plt.ylim(data['y_train'].min() - margin, data['y_train'].max() + margin)
plt.legend(loc='upper right', prop={'size': 20})
plt.title('Fitted polynomial regressions with L2 reg')
plt.xlabel('x')
plt.ylabel('y')
plt.show() 

В результате даже 13-ая степень ведет себя так, как мы ожидаем. Графики немного сгладились, хотя мы все равно наблюдаем небольшое переобучение на степенях выше третьей, что выражается в интерполяции данных в правой части графика.

for ix, t in enumerate(w_list_l2):
    degree, w = t
    fig, ax = plt.subplots()
    plt.bar(range(max(degree_list) + 1), np.hstack((w, [0]*(max(degree_list) - w.shape[0] + 1))), color=colors[ix])
    plt.title('Magnitude of fitted LR parameters for poly:%i with L2 reg' % degree)
    plt.xlabel('degree')
    plt.ylabel('value of w')
    ax.set_xticks(np.array(range(max(degree_list) + 1)) + 0.5)
    ax.set_xticklabels(range(max(degree_list) + 1))
    plt.show()

Амплитуда коэффициентов также изменилась, хотя скакать в разные стороны они не перестали. Мы помним, что полином третьей степени должен лучше всего описывать наши данные, хотелось бы, чтобы в результате регуляризации все коэффициенты при полиномиальных признаках степени выше третьей были равны нулю. И, оказывается, есть и такой регуляризатор.

регуляризация

lmbd = 1
degree = 13
dlist = [[1]*data['x_train'].shape[0]] + map(lambda n: data['x_train']**n, range(1, degree + 1))
X = np.array(dlist).T

# функция для вычисления среднеквадратичное ошибки
def mse(u, v):
    return ((u - v)**2).sum()/u.shape[0]

# начальное приближение
w = np.array([-1.0] * X.shape[1])
# максимальное количество итераций
n_iter = 20
# сделаем скорость обучения очень маленькой, на всякий случай 
lr = 0.00000001
loss = []
for ix in range(n_iter):
    w -= lr*(np.dot(np.dot(X, w) - data['y_train'], X)/X.shape[0] + lmbd*np.sign(w))
    y_hat = np.dot(X, w)
    loss.append(mse(data['y_train'], y_hat))
    print loss[-1]

Даже при такой небольшой скорости обучения, ошибка все равно растет и очень даже стремительно. Причина в том, что каждый признак измеряется в разных масштабах, от небольших чисел у полиномиальных признаков 1–2 степени, до огромных при 12–13 степени. Для того чтобы итеративный процесс сошелся, необходимо либо выбрать экстремально мелкую скорость обучения, либо каким-то образом нормализовать признаки. Применим следующее преобразование к признакам и попробуем запустить процесс еще раз:

Такое преобразование называется стандартизацией, распределение каждого признака теперь имеет нулевое матожидание и единичную дисперсию.

lmbd = 1
degree = 13
dlist = [[1]*data['x_train'].shape[0]] + map(lambda n: data['x_train']**n, range(1, degree + 1))
X = np.array(dlist).T
# вычислим выборочное среднее каждого признака
x_mean = X.mean(axis=0)
# вычислим выборочное стандартное отклонение признаков
x_std = X.std(axis=0)
# применим преобразование
X = (X - x_mean)/x_std
X[:, 0] = 1.0

w = np.array([-1.0] * X.shape[1])
n_iter = 100
lr = 0.1
loss = []
for ix in range(n_iter):
    w -= lr*(np.dot(np.dot(X, w) - data['y_train'], X)/X.shape[0] + lmbd*np.sign(w))
    y_hat = np.dot(X, w)
    loss.append(mse(data['y_train'], y_hat))

plt.plot(loss)
plt.title('Train error')
plt.xlabel('Iteration')
plt.ylabel('MSE')
plt.show()

Все стало сильно лучше.

Нарисуем теперь все графики:

degree_list = [1, 2, 3, 10, 12, 13]
cmap = plt.get_cmap('jet')
colors = [cmap(i) for i in np.linspace(0, 1, len(degree_list))]

margin = 0.3
plt.plot(data['support'], data['values'], 'b--', alpha=0.5, label='manifold')
plt.scatter(data['x_train'], data['y_train'], 40, 'g', 'o', alpha=0.8, label='data')

def mse(u, v):
    return ((u - v)**2).sum()/u.shape[0]

def fit_lr_l1(X, y, lmbd, n_iter=100, lr=0.1):
    w = np.array([-1.0] * X.shape[1])
    loss = []
    for ix_iter in range(n_iter):
        w -= lr*(np.dot(np.dot(X, w) - y, X)/X.shape[0] +lmbd*np.sign(w))
        y_hat = np.dot(X, w)
        loss.append(mse(y, y_hat))
    return w, y_hat, loss
    
w_list_l1 = []
for ix, degree in enumerate(degree_list):
    dlist = [[1]*data['x_train'].shape[0]] + map(lambda n: data['x_train']**n, range(1, degree + 1))
    X = np.array(dlist).T
    x_mean = X.mean(axis=0)
    x_std = X.std(axis=0)
    X = (X - x_mean)/x_std
    X[:, 0] = 1.0
    w, y_hat, loss = fit_lr_l1(X, data['y_train'], lmbd=0.05)
    w_list_l1.append((degree, w))
    plt.plot(data['x_train'], y_hat, color=colors[ix], label='poly degree: %i' % degree)

plt.xlim(data['x_train'].min() - margin, data['x_train'].max() + margin)
plt.ylim(data['y_train'].min() - margin, data['y_train'].max() + margin)
plt.legend(loc='upper right', prop={'size': 20})
plt.title('Fitted polynomial regressions with L1 reg')
plt.xlabel('x')
plt.ylabel('y')
plt.show() 

Если посмотреть на коэффициенты, мы увидим, что большая часть из них близка к нулю (то, что у 13-ой степени коэффициент совсем не нулевой, можно списать на шум и малое количество примеров в обучающей выборке; так же стоит помнить, что теперь все признаки измеряются в одинаковых шкалах).

for ix, t in enumerate(w_list_l1):
    degree, w = t
    fig, ax = plt.subplots()
    plt.bar(range(max(degree_list) + 1), np.hstack((w, [0]*(max(degree_list) - w.shape[0] + 1))), color=colors[ix])
    plt.title('Magnitude of fitted LR parameters for poly:%i with L1 reg' % degree)
    plt.xlabel('degree')
    plt.ylabel('value of w')
    ax.set_xticks(np.array(range(max(degree_list) + 1)) + 0.5)
    ax.set_xticklabels(range(max(degree_list) + 1))
    plt.show()

Описанный способ построения регрессии называется LASSO регрессия. Очень хотелось бы думать, что дядька на коне бросает веревку и ворует коэффициенты, а на их месте остается нуль. Но нет, LASSO = least absolute shrinkage and selection operator.

Байесовская интерпретация линейной регрессии

•  — априорные ожидания (prior): насколько правдоподобна гипотеза перед наблюдением данных;
•  — правдоподобие (likelihood): насколько правдоподобны данные при условии того, что гипотеза верна;
•  — маргинальная вероятность (marginal probability или evidence): вероятность данных, усредненная по всевозможным гипотезам;
•  — апостериорное распределение (posterior): насколько правдоподобна гипотеза при наблюдаемых данных.

В статистике обычно ищут точечную оценку максимума правдоподобия (ML = maximum likelihood):

В то время как в байесовом подходе интересуются апостериорным распределением:

Часто получается так, что интеграл, полученный в результате байесового вывода, крайне нетривиален (в случае линейной регрессии это, к счастью, не так), и тогда нужна точечная оценка. Тогда мы интересуемся максимумом апостериорного распределения (MAP = maximum a posteriori):

Давайте сравним ML и MAP гипотезы для линейной регрессии, это даст нам четкое понимание смысла регуляризаций. Будем считать, что все объекты из обучающей выборки были взяты из общей популяции независимо и равномерно распределенно. Это позволит нам записать совместную вероятность данных (правдоподобие) в виде:

А также будем считать, что целевая переменная подчиняется следующему закону:

data = generate_wave_set(1000, 100)
X = np.vstack((np.ones(data['x_train'].shape[0]), data['x_train'])).T
w = np.dot(np.dot(np.linalg.inv(np.dot(X.T, X)), X.T), data['y_train'])

w0_support = np.linspace(-3, 3, 1000)
w1_support = np.linspace(-3, 3, 1000)
# create cartesian product of parameters
wx_space = list(it.product(w0_support, w1_support))
w0, w1 = zip(*wx_space)
# calculate MSE on dataset for each pairs of parameters
y = ((data['y_train'][:, np.newaxis] - np.dot(X, np.array(wx_space).T))**2).mean(axis=0)

plt.hexbin(w0, w1, C=y**(0.2), cmap=cm.jet_r, bins=None)
plt.axvline(0, color='black', linestyle='-', label='origin')
plt.axhline(0, color='black', linestyle='-')
plt.axvline(w[0], color='w', linestyle='--', label='ML solution')
plt.axhline(w[1], color='w', linestyle='--')
plt.axes().set_aspect('equal', 'datalim')
plt.title('ML solution')
plt.xlabel('w_0')
plt.ylabel('w_1')
plt.legend(loc='upper left', prop={'size': 20})
plt.show()

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

Найдем MAP оценку параметров линейной регрессии, для этого придется задать какое-нибудь априорное распределение на параметры модели. Пусть для начала это будет опять нормальное распределение: .

x = np.linspace(-5, 5, 1000)
for scale in np.linspace(0.5, 1.4, 7):
    plt.plot(x, norm.pdf(x, scale=scale), label='scale=%0.2f' % scale)
    
plt.legend(loc='upper right', prop={'size': 20})
plt.title('Normal distribution with different scale parameter')
plt.show()

Тогда апостериорное распределение примет вид:

© Habrahabr.ru