Перевод. Оригинал: Introduction to Regression Splines (with Python codes)
Введение
Будучи новичком в мире наук о данных, первым алгоритмом, который я изучил, была линейная регрессия. Я применил его к различным наборам данных и оценил как его преимущества, так и ограничения.
Предполагалось, что между зависимой и независимой переменными существует линейная связь, что в действительности редко имело место. В качестве улучшения по сравнению с этой моделью я попробовал полиномиальную регрессию, которая давала лучшие результаты (большую часть времени). Но использование полиномиальной регрессии для наборов данных с высокой степенью изменчивости может привести к переобучению.
Моя модель всегда становилась слишком гибкой, и плохо работала с другими данными. Затем я наткнулся на другой нелинейный подход, известный как сплайновая регрессия. Она использует комбинацию линейных/полиномиальных функций, чтобы аппроксимировать данные.
В этой статье мы рассмотрим некоторые основы линейной и полиномиальной регрессии и подробно изучим значение сплайнов и их реализацию в Python.
Примечание. Для полного понимания концепций, рассматриваемых в этой статье, необходимы знания о линейной и полиномиальной регрессии. Вы можете узнать больше о них здесь.
Понимание данных
Чтобы понять концепции, мы будем работать с набором данных для прогнозирования заработной платы, который вы можете скачать здесь (он был взят из популярной книги «Introduction to Statistical learning»).
Наш набор данных содержит такую информацию, как ID, год, возраст, пол, семейное положение, раса, образование, регион, тип работы, здоровье, медицинское страхование, журнал заработной платы и заработное платы различных сотрудников. Чтобы сфокусироваться на сплайн-регрессии, я буду использовать в качестве независимой переменной только «возраст» для прогнозирования заработной платы (зависимая переменная).
Давайте начнем работать с данными.
# импорт модулей import pandas as pd import numpy as np import statsmodels.api as sm import matplotlib.pyplot as plt %matplotlib inline # читаем набор данных data = pd.read_csv("Wage.csv") data.head()
data_x = data['age'] data_y = data['wage'] # делим набор данных на обучающую и тестовую выборки from sklearn.model_selection import train_test_split train_x, valid_x, train_y, valid_y = train_test_split(data_x, data_y, test_size=0.33, random_state = 1) # Визуализируем отношения возраста и заработной платы import matplotlib.pyplot as plt plt.scatter(train_x, train_y, facecolor='None', edgecolor='k', alpha=0.3) plt.show()
Что вы думаете о приведенном графике рассеяния? Это положительная, отрицательная или отсутствие корреляции?
Введение в линейную регрессию
Линейная регрессия является самым простым и наиболее широко используемым статистическим методом для прогнозного моделирования. Это алгоритм с учителем для решения регрессионных задач.
Она называется линейной, поскольку устанавливает линейную связь между зависимой и независимой переменными. Это дает нам линейное уравнение, подобное приведенному ниже, где мы имеем наши признаки в качестве независимых переменных с коэффициентами:
Здесь Y - наша зависимая переменная, X - независимые переменные, а бета - коэффициенты. Коэффициенты - это веса, присвоенные элементам. Они означают важность каждого из признаков. Например, если результат уравнения сильно зависит от одного признака (X1) по сравнению с любым другим признаком, это означает, что коэффициент/вес признака (X1) будет иметь более высокую величину по сравнению с любым другим признаком.
Итак, давайте попробуем понять линейную регрессию с одним признаком, то есть только с одной независимой переменной. Это называется простая линейная регрессия. Следовательно, наше уравнение становится следующим:
Поскольку мы используем только «возраст» для прогнозирования «заработной платы» сотрудников, мы будем реализовывать простую линейную регрессию на обучающей выборке и вычислять погрешность (RMSE) по тестовой выборке.
from sklearn.linear_model import LinearRegression # обучаем модель регрессии x = train_x.reshape(-1,1) model = LinearRegression() model.fit(x,train_y) print(model.coef_) print(model.intercept_)
-> array([0.72190831])
-> 80.65287740759283
# Прогноз на тестовой выборке valid_x = valid_x.reshape(-1,1) pred = model.predict(valid_x) # Визуализация # Мы будем использовать 70 точек valid_x xp = np.linspace(valid_x.min(),valid_x.max(),70) xp = xp.reshape(-1,1) pred_plot = model.predict(xp) plt.scatter(valid_x, valid_y, facecolor='None', edgecolor='k', alpha=0.3) plt.plot(xp, pred_plot) plt.show()
Теперь мы можем рассчитать RMSE по прогнозам.
from sklearn.metrics import mean_squared_error from math import sqrt rms = sqrt(mean_squared_error(valid_y, pred)) print(rms) -> 40.436
Из вышеприведенного графика можно сделать вывод, что линейная регрессия не охватывает все доступные сигналы и не является лучшим методом для получения этого прогноза заработной платы.
Хотя линейные модели относительно просты в описании и реализации и имеют преимущества перед другими подходами с точки зрения интерпретации, они имеют существенные ограничения с точки зрения предсказательной силы, потому что они предполагают линейную связь между зависимой и независимой переменными, которая почти всегда является приблизительной, а иногда и не является таковой.
В других методах, которые мы увидим ниже, мы отбросим предположение о линейности, все еще пытаясь поддерживать как можно большую интерпретируемость. Мы сделаем это путем изучения очень простых расширений линейных моделей, таких как полиномиальная регрессия и ступенчатые функции, а также более сложных подходов, таких как сплайны.
Улучшение по сравнению с линейной регрессией: полиномиальная регрессия
Рассмотрим эти визуализации:
Графики выше, кажется, используют намного больше сигналов между заработной платой и возрастом по сравнению с линейным графиком. Эти графики не являются линейными по форме, следовательно, они используют нелинейное уравнение вместо линейного уравнения для установления взаимосвязи между возрастом и заработной платой. Этот тип регрессии, который использует нелинейные функции, называется полиномиальной регрессией.
Полиномиальная регрессия расширяет линейную модель, добавляя дополнительные предикторы, получаемые путем возведения в степень каждого из исходных предикторов. Например, кубическая регрессия использует три переменные:
в качестве предикторов. Этот подход обеспечивает простой способ обеспечить нелинейное соответствие данных.
в качестве предикторов. Этот подход обеспечивает простой способ обеспечить нелинейное соответствие данных.
Стандартный метод расширения линейной регрессии до нелинейной связи между зависимой и независимой переменными состоял в том, чтобы заменить линейную модель полиномиальной функцией.
Когда мы увеличиваем значение степени, полученная кривая содержит большие колебания, которые приведут к чрезмерно гибким формам графика. Такие кривые приводят к переобучению.
# Генерируем веса для полиномиальной функции степени 2 weights = np.polyfit(train_x, train_y, 2) print(weights) -> array([ -0.05194765, 5.22868974, -10.03406116]) # Генерируем модель с данными весами model = np.poly1d(weights) # Прогноз на тестовой выборке pred = model(valid_x) # Строим график только для 70 наблюдений xp = np.linspace(valid_x.min(),valid_x.max(),70) pred_plot = model(xp) plt.scatter(valid_x, valid_y, facecolor='None', edgecolor='k', alpha=0.3) plt.plot(xp, pred_plot) plt.show()
Точно так же мы можем построить полиномиальные кривые для разных значений степени.
К сожалению, полиномиальная регрессия также имеет немало проблем. По мере того, как мы увеличиваем сложность формулы, число признаков также увеличивается, что иногда трудно рассчитать. Кроме того, полиномиальная регрессия имеет тенденцию к резкому переобучению даже на этом простом одномерном наборе данных.
Есть и другие проблемы с полиномиальной регрессией. Например, она по своей природе нелокальна, то есть изменение значения одной точки в обучающем наборе может повлиять на подгонку полинома для точек данных, которые находятся очень далеко. Следовательно, чтобы избежать использования полинома высокой степени для всего набора данных, мы можем заменить его несколькими различными полиномиальными функциями малой степени.
Сплайновая регрессия и ее реализации
Чтобы преодолеть недостатки полиномиальной регрессии, мы можем использовать улучшенный метод регрессии, который вместо построения одной модели для всего набора данных делит набор данных на несколько блоков и подгоняет каждый блок отдельной моделью. Такая техника называется сплайн-регрессией.
Сплайн-регрессия - один из наиболее важных методов нелинейной регрессии. При полиномиальной регрессии мы генерировали новые функции, используя различные полиномиальные функции для существующих признаков, которые навязывали глобальную структуру для набора данных. Чтобы преодолеть это, мы можем разделить набор данных на отдельные части и подогнать линейные или полиномиальные функции низкой степени на каждой из этих частей.
Точки, где происходит разделение, называются узлами. Функции, которые мы можем использовать для моделирования каждого куска/бина, называются кусочными функциями. Существуют различные кусочные функции, которые мы можем использовать для подгонки этих отдельных участков.
В следующих нескольких подразделах мы узнаем о некоторых из этих кусочных функций.
Кусочно-шаговые функции
Одной из наиболее распространенных кусочных функций является шаговая функция. Шаговая функция - это функция, которая остается постоянной в пределах интервала. Мы можем подогнать отдельные пошаговые функции к каждой из разделенных частей, чтобы избежать навязывания глобальной структуры. Здесь мы разбиваем диапазон X на интервалы и устанавливаем разные константы для каждого интервала.
То есть, мы создаем точки разрезания диапазона X C1, C2,. , , , Ck, а затем строим K+1 новых переменных.
где I() - это индикаторная функция, которая возвращает 1, если условие истинно, и 0 в противном случае. Например, I(cK ≤ X) равно 1, если cK ≤ X, в противном случае оно равно 0. Для данного значения X не более одного из C1, C2,. , ., CK может быть ненулевым, так как X должен лежать только в одном из интервалов.
# Делим данные на 4 интервала df_cut, bins = pd.cut(train_x, 4, retbins=True, right=True) df_cut.value_counts(sort=False) -> (17.938, 33.5] 504 (33.5, 49.0] 941 (49.0, 64.5] 511 (64.5, 80.0] 54 Name: age, dtype: int64
df_steps = pd.concat([train_x, df_cut, train_y], keys=['age','age_cuts','wage'], axis=1)
# Создаем фиктивные переменные для возрастных групп
df_steps_dummies = pd.get_dummies(df_cut)
df_steps_dummies.head()
df_steps_dummies.columns = ['17.938-33.5','33.5-49','49-64.5','64.5-80'] # Подгоняем обобщенные линейные модели fit3 = sm.GLM(df_steps.wage, df_steps_dummies).fit() # Разделяем тестовую выборку на 4 интервала bin_mapping = np.digitize(valid_x, bins) X_valid = pd.get_dummies(bin_mapping) # Удаляем выбросы X_valid = pd.get_dummies(bin_mapping).drop([5], axis=1) # Прогноз pred2 = fit3.predict(X_valid) # Рассчитываем RMSE from sklearn.metrics import mean_squared_error from math import sqrt rms = sqrt(mean_squared_error(valid_y, pred2)) print(rms) ->39.9 # Строим график только для 70 наблюдений xp = np.linspace(valid_x.min(),valid_x.max()-1,70) bin_mapping = np.digitize(xp, bins) X_valid_2 = pd.get_dummies(bin_mapping) pred2 = fit3.predict(X_valid_2)
# Визуализация
fig, (ax1) = plt.subplots(1,1, figsize=(12,5))
fig.suptitle('Piecewise Constant', fontsize=14)
# График рассеяния с кривой полиномиальной регрессии
ax1.scatter(train_x, train_y, facecolor='None', edgecolor='k', alpha=0.3)
ax1.plot(xp, pred2, c='b')
ax1.set_xlabel('age')
ax1.set_ylabel('wage')
plt.show()
Разделение на интервалы имеет свои очевидные концептуальные проблемы. Наиболее важно, что мы ожидаем, что большинство изучаемых нами явлений будет постоянно зависеть от входных данных. Интервальная регрессия не создает непрерывных функций предиктора, поэтому в большинстве случаев мы не можем ожидать никакой связи между входом и выходом.
Например, на приведенном выше графике мы видим, что на первом интервале явно отсутствует тенденция к росту заработной платы с возрастом.
Базисные функции
Чтобы отразить нелинейность в регрессионных моделях, нам нужно преобразовать некоторые или все предикторы. Чтобы избежать необходимости рассматривать каждый предиктор как линейный, мы хотим применить к ним очень общее семейство преобразований. Семейство должно быть достаточно гибким, чтобы адаптироваться (при подгонке моделей) к широкому спектру форм, но не должно быть слишком гибким, чтобы не допустить переобучения.
Эта концепция семейства преобразований, которые могут объединяться для захвата общих форм, называется базисной функцией. В этом случае нашими объектами являются функции: b1(X), b2(X),. , , , BK(X).
Вместо подгонки линейной модели по X мы подгоняем следующую модель:
Теперь мы рассмотрим очень распространенный выбор для базисной функции: кусочно-полиномиальные функции.
Кусочно-полиномиальные функции
Вместо того, чтобы подгонять постоянную функцию по разным интервалам в диапазоне X, кусочно-полиномиальная регрессия включает в себя подгонку отдельных полиномов низкой степени к различным интервалам X. Поскольку мы используем более низкие степени полиномов, мы не наблюдаем больших колебаний кривой вокруг данных.
Например, кусочно-квадратичный полином работает, подгоняя квадратное уравнение регрессии:
где коэффициенты β0, β1 и β2 различаются в разных интервалах диапазона X. Кусочно-кубический многочлен с одним узлом в точке c принимает следующую форму:
Другими словами, мы подгоняем к данным две разные полиномиальные функции: одну на подмножестве наблюдений с xi
Первая полиномиальная функция имеет коэффициенты β01, β11, β21, β31, а вторая имеет коэффициенты β02, β12, β22, β32. Каждая из этих полиномиальных функций может быть подобрана с использованием метода наименьших квадратов.
Помните, что это семейство полиномиальных функций имеет 8 степеней свободы, 4 для каждого полинома (так как есть 4 переменные).
Использование большего количества узлов приводит к более гибкому кусочному полиному, так как мы используем разные функции для каждого интервала. Эти функции зависят только от распределения данных этого конкретного интервала. В общем, если мы поместим K различных узлов во всем диапазоне X, мы получим в итоге K + 1 различных кубических многочленов. Мы можем использовать любой многочлен низкой степени, чтобы соответствовать этим отдельным интервалам. Или, вместо этого, мы можем использовать кусочно-линейные функции. Фактически использованные выше шаговые функции на самом деле являются кусочными полиномами степени 0.
Теперь мы рассмотрим некоторые необходимые условия и ограничения, которые должны соблюдаться при формировании кусочных полиномов.
Ограничения и сплайны
Мы должны быть осторожны при использовании кусочных полиномов, поскольку существуют различные ограничения, которым мы должны следовать. Посмотрите на изображение ниже:
Мы можем столкнуться с определенными ситуациями, когда многочлены не являются непрерывными в узле. Такого условия следует избегать, потому что семейство полиномов в целом должно генерировать уникальный вывод для каждого входа.
Из приведенного выше изображения видно, что на первом узле выводятся два разных значения. Таким образом, чтобы избежать этого, мы должны добавить дополнительное ограничение/условие, чтобы полиномы по обе стороны от узла были непрерывными в узле.
Теперь после добавления этого ограничения мы получаем непрерывное семейство полиномов. Но выглядит ли это идеально? Прежде чем читать дальше, подумайте, чего здесь не хватает.
Похоже, гладкость у узлов пока отсутствует. Таким образом, чтобы сгладить многочлены в узлах, мы добавляем дополнительное ограничение/условие: первая производная обоих многочленов должна быть одинаковой. Следует отметить одну вещь: каждое ограничение, которое мы накладываем на кусочно-кубические полиномы, эффективно освобождает одну степень свободы, поскольку мы уменьшаем сложность получающейся кусочно-полиномиальной модели. Поэтому на приведенном выше графике мы используем только 10 степеней свободы вместо 12.
После наложения ограничения о равенстве первой производной, получаем приведенный выше график. Этот график использует 8 степеней свободы вместо 12, поскольку налагаются два ограничения. Хотя приведенный выше сюжет выглядит лучше, все же есть еще некоторые возможности для улучшения. Теперь мы наложим дополнительное ограничение: вторые производные обоих многочленов в узле должны быть одинаковыми.
Этот график кажется идеальным для нашего исследования. Он использует 6 степеней свободы вместо 12. Такой кусочный многочлен степени m с m-1 непрерывными производными называется сплайном. Следовательно, мы построили кубический сплайн на графике выше. Мы можем построить сплайн любой степени с m-1 непрерывными производными.
Кубические и натуральные кубические сплайны
Кубический сплайн является кусочным полиномом с набором дополнительных ограничений (непрерывность, непрерывность первой производной и непрерывность второй производной). В общем, кубический сплайн с K узлами использует кубический сплайн с 4 + K степенями свободы. Редко есть какая-либо веская причина выйти за пределы кубических сплайнов (если только не интересоваться гладкими производными).
from patsy import dmatrix import statsmodels.api as sm import statsmodels.formula.api as smf # Генерация кубического сплайна transformed_x = dmatrix("bs(train, knots=(25,40,60), degree=3, include_intercept=False)", {"train": train_x},return_type='dataframe') # Подгонка обобщенной линейной модели на преобразованном наборе данных fit1 = sm.GLM(train_y, transformed_x).fit() # Генерируем кубический сплайн с 4 узлами transformed_x2 = dmatrix("bs(train, knots=(25,40,50,65),degree =3, include_intercept=False)", {"train": train_x}, return_type='dataframe') # Подгонка обобщенной линейной модели на преобразованном наборе данных fit2 = sm.GLM(train_y, transformed_x2).fit() # Прогнозы на обоих сплайнах pred1 = fit1.predict(dmatrix("bs(valid, knots=(25,40,60), include_intercept=False)", {"valid": valid_x}, return_type='dataframe')) pred2 = fit2.predict(dmatrix("bs(valid, knots=(25,40,50,65),degree =3, include_intercept=False)", {"valid": valid_x}, return_type='dataframe')) # Рассчитываем значения RMSE rms1 = sqrt(mean_squared_error(valid_y, pred1)) print(rms1) -> 39.4 rms2 = sqrt(mean_squared_error(valid_y, pred2)) print(rms2) -> 39.3 # Строим график только для 70 наблюдений xp = np.linspace(valid_x.min(),valid_x.max(),70) # Делаем прогнозы pred1 = fit1.predict(dmatrix("bs(xp, knots=(25,40,60), include_intercept=False)", {"xp": xp}, return_type='dataframe')) pred2 = fit2.predict(dmatrix("bs(xp, knots=(25,40,50,65),degree =3, include_intercept=False)", {"xp": xp}, return_type='dataframe')) # Строим сплайны и каналы ошибок plt.scatter(data.age, data.wage, facecolor='None', edgecolor='k', alpha=0.1) plt.plot(xp, pred1, label='Specifying degree =3 with 3 knots') plt.plot(xp, pred2, color='r', label='Specifying degree =3 with 4 knots') plt.legend() plt.xlim(15,85) plt.ylim(0,350) plt.xlabel('age') plt.ylabel('wage') plt.show()
Мы знаем, что поведение многочленов, которые соответствуют данным, имеет тенденцию быть неустойчивым вблизи границ. Такая изменчивость может быть опасной. Эти проблемы также относятся и к сплайнам. Полиномы, выходящие за пределы граничных узлов, ведут себя еще более дико, чем соответствующие глобальные полиномы в этой области. Чтобы сгладить полином за граничными узлами, мы будем использовать специальный тип сплайна, известный как естественный сплайн (Natural Spline).
Естественный кубический сплайн добавляет дополнительные ограничения, а именно, что функция является линейной за граничными узлами. Это ограничивает кубическую и квадратичную части там нулем, что уменьшает количество степеней свободы на 2 для каждого случая. Это 2 степени свободы на каждом из двух концов кривой, уменьшая K + 4 до K.
# Генерируем естественный кубический сплайн transformed_x3 = dmatrix("cr(train,df = 3)", {"train": train_x}, return_type='dataframe') fit3 = sm.GLM(train_y, transformed_x3).fit() # Прогноз на проверочной выборке pred3 = fit3.predict(dmatrix("cr(valid, df=3)", {"valid": valid_x}, return_type='dataframe')) # Рассчитываем значение RMSE rms = sqrt(mean_squared_error(valid_y, pred3)) print(rms) -> 39.44 # Строим график только для 70 наблюдений xp = np.linspace(valid_x.min(),valid_x.max(),70) pred3 = fit3.predict(dmatrix("cr(xp, df=3)", {"xp": xp}, return_type='dataframe')) # Строим график сплайна plt.scatter(data.age, data.wage, facecolor='None', edgecolor='k', alpha=0.1) plt.plot(xp, pred3,color='g', label='Natural spline') plt.legend() plt.xlim(15,85) plt.ylim(0,350) plt.xlabel('age') plt.ylabel('wage') plt.show()
Выбор количества и местоположения узлов
Когда мы подгоняем сплайн, где мы должны разместить узлы? Одним из потенциальных мест может быть область высокой изменчивости, потому что в этих областях коэффициенты полинома могут быстро меняться. Следовательно, одним из вариантов является размещение большего количества узлов в тех местах, где, по нашему мнению, функция может изменяться наиболее быстро, и размещение меньшего количества узлов в тех местах, где она кажется более стабильной.
Хотя эта опция может работать хорошо, на практике обычно принято размещать узлы равномерно. Один из способов сделать это - указать желаемые степени свободы, а затем заставить программное обеспечение автоматически размещать соответствующее количество узлов в одинаковых квантилях данных.
Другой вариант - опробовать разное количество узлов и посмотреть, какой вариант дает наилучшую кривую.
Более объективный подход заключается в использовании кросс-валидации. С помощью этого метода: мы удаляем часть данных, подгоняем сплайн с определенным количеством узлов к оставшимся данным, а затем используем сплайн, чтобы сделать прогноз для удаленной части данных.
Мы повторяем этот процесс несколько раз, пока каждое наблюдение не будет пропущено хотя один раз, а затем вычисляем общее перекрестное значение RMSE. Эту процедуру можно повторить для разного количества K узлов. Затем выбирается значение K, дающее наименьшее RMSE.
Сравнение сплайновой регрессии с полиномиальной регрессией
Сплайновая регрессия часто дает лучшие результаты, чем полиномиальная регрессия. Это связано с тем, что в отличие от полиномов высоких степеней, которые необходимо использовать для получения гибкой подгонки, сплайны обеспечивают гибкость, увеличивая количество узлов, но сохраняя фиксированную степень.
Как правило, этот подход дает более стабильные оценки. Сплайны также позволяют нам размещать больше узлов и, следовательно, дают большую гибкость в областях, где функция быстро меняется, и меньше узлов, где функция кажется более стабильной. Дополнительная гибкость полинома приводит к нежелательным результатам на границах, тогда как естественный кубический сплайн все же обеспечивает разумное соответствие данным.
Комментариев нет:
Отправить комментарий