1.1. Линейные модели

1.1.1. Метод наименьших квадратов

LinearRegression подбирает линейную модель с коэффициентами \(w = (w_1, ..., w_p)\), минимизируя сумму квадратов остатков между наблюдаемыми целевыми значениями в наборе данных и значениями, предсказанными линейной аппроксимацией. Математически решается задача следующего вида:

LinearRegression принимает в метод fit аргументы X, y, sample_weight и сохраняет коэффициенты \(w\) линейной модели в атрибутах coef_ и intercept_:

>>> from sklearn import linear_model
>>> reg = linear_model.LinearRegression()
>>> reg.fit([[0, 0], [1, 1], [2, 2]], [0, 1, 2])
LinearRegression()
>>> reg.coef_
array([0.5, 0.5])
>>> reg.intercept_
0.0

Оценки коэффициентов методом наименьших квадратов основываются на независимости признаков. Если признаки коррелируют и некоторые столбцы матрицы плана \(X\) приблизительно линейно зависимы, матрица плана становится почти вырожденной. В результате оценка методом наименьших квадратов становится крайне чувствительной к случайным ошибкам в наблюдаемой целевой переменной, что приводит к большой дисперсии. Такая ситуация, называемая мультиколлинеарностью, может возникнуть, например, если данные собирались без экспериментального плана.

Примеры

• Метод наименьших квадратов и регрессия Ridge

1.1.2. Регрессия и классификация Ridge

1.1.2.1. Регрессия

Регрессия Ridge решает некоторые проблемы метода наименьших квадратов, накладывая штраф на величину коэффициентов. Коэффициенты Ridge минимизируют сумму квадратов остатков со штрафом:

\[\min_{w} || X w - y||_2^2 + \alpha ||w||_2^2\]

Параметр сложности \(\alpha \geq 0\) определяет степень сжатия: чем больше значение \(\alpha\), тем сильнее сжатие, а значит, коэффициенты становятся устойчивее к коллинеарности.

Как и другие линейные модели, Ridge принимает в метод fit массивы X, y и сохраняет коэффициенты \(w\) линейной модели в атрибуте coef_:

>>> from sklearn import linear_model
>>> reg = linear_model.Ridge(alpha=.5)
>>> reg.fit([[0, 0], [0, 0], [1, 1]], [0, .1, 1])
Ridge(alpha=0.5)
>>> reg.coef_
array([0.34545455, 0.34545455])
>>> reg.intercept_
np.float64(0.13636)

Обратите внимание, что класс Ridge позволяет пользователю указать, что решатель должен выбираться автоматически, задав solver="auto". Если этот параметр указан, Ridge выберет один из решателей: "lbfgs", "cholesky" или "sparse_cg". Ridge будет проверять условия из следующей таблицы сверху вниз. Если условие выполняется, выбирается соответствующий решатель.

Решатель	Условие
‘lbfgs’	Указан параметр positive=True.
‘cholesky’	Входной массив X не является разреженным.
‘sparse_cg’	Ни одно из перечисленных выше условий не выполняется.

Примеры

• Метод наименьших квадратов и регрессия Ridge
• График зависимости коэффициентов Ridge от регуляризации
• Распространённые ошибки при интерпретации коэффициентов линейных моделей
• Зависимость коэффициентов Ridge от L2-регуляризации

>>> import numpy as np
>>> from sklearn import linear_model
>>> reg = linear_model.RidgeCV(alphas=np.logspace(-6, 6, 13))
>>> reg.fit([[0, 0], [0, 0.1], [1, 1]], [0, -0.1, 1])
RidgeCV(alphas=array([1.e-06, 1.e-05, 1.e-04, 1.e-03, 1.e-02, 1.e-01, 1.e+00, 1.e+01,
      1.e+02, 1.e+03, 1.e+04, 1.e+05, 1.e+06]))
>>> reg.alpha_
np.float64(0.1)

1.1.3. Лассо

Lasso — это линейная модель, которая оценивает разреженные коэффициенты, то есть способна точно обнулять коэффициенты. В некоторых случаях она полезна благодаря склонности предпочитать решения с меньшим количеством ненулевых коэффициентов, фактически уменьшая число признаков, от которых зависит данное решение. По этой причине Лассо и его варианты играют фундаментальную роль в области сжатого зондирования. При определённых условиях метод способен восстановить точный набор ненулевых коэффициентов (см. Сжатое зондирование: реконструкция томографии с априорным распределением L1 (Лассо)).

Математически это линейная модель с добавленным членом регуляризации. Минимизируемая целевая функция имеет вид:

Таким образом, оценка Лассо решает задачу наименьших квадратов с добавленным штрафом \(\alpha ||w||_1\), где \(\alpha\) — константа, а \(||w||_1\) — норма \(\ell_1\) вектора коэффициентов.

В реализации класса Lasso для подгонки коэффициентов используется алгоритм покоординатного спуска. Другую реализацию см. в разделе Регрессия наименьшего угла:

>>> from sklearn import linear_model
>>> reg = linear_model.Lasso(alpha=0.1)
>>> reg.fit([[0, 0], [1, 1]], [0, 1])
Lasso(alpha=0.1)
>>> reg.predict([[1, 1]])
array([0.8])

Функция lasso_path полезна для задач низкого уровня, поскольку вычисляет коэффициенты вдоль всего пути возможных значений.

Примеры

• Модели на основе L1 для разреженных сигналов
• Сжатое зондирование: реконструкция томографии с априорным распределением L1 (Лассо)
• Распространённые ошибки при интерпретации коэффициентов линейных моделей
• Выбор модели Лассо: AIC-BIC / кросс-валидация

1.1.3.2. Настройка параметра регуляризации

Параметр alpha определяет степень разреженности оценённых коэффициентов.

1.1.3.2.1. Использование кросс-валидации

scikit-learn предоставляет объекты для настройки параметра alpha метода Лассо с помощью кросс-валидации: LassoCV и LassoLarsCV. LassoLarsCV основан на алгоритме регрессии наименьшего угла, описанном ниже.

Для наборов данных высокой размерности с большим количеством коллинеарных признаков чаще всего предпочтительнее использовать LassoCV. Однако LassoLarsCV позволяет исследовать более подходящие значения параметра alpha и часто работает быстрее, чем LassoCV, если число образцов очень мало по сравнению с числом признаков.

1.1.3.2.2. Выбор модели на основе информационных критериев

В качестве альтернативы оцениватель LassoLarsIC предлагает использовать информационный критерий Акаике (AIC) и байесовский информационный критерий (BIC). Это более дешёвая с вычислительной точки зрения альтернатива для поиска оптимального значения alpha, поскольку путь регуляризации вычисляется только один раз, а не k+1 раз, как при использовании k-блочной кросс-валидации.

Эти критерии вычисляются на обучающем наборе, используемом для подгонки модели. Вкратце, они штрафуют чрезмерно оптимистичные оценки различных моделей Лассо с учётом их гибкости (см. раздел «Математические подробности» ниже).

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

Примеры

• Выбор модели Лассо: AIC-BIC / кросс-валидация
• Выбор модели Лассо с помощью информационных критериев

1.1.3.2.3. Критерии AIC и BIC

Определение AIC (и, следовательно, BIC) в литературе может различаться. В этом разделе приводятся дополнительные сведения о критерии, вычисляемом в scikit-learn.

Математические подробности

Критерий AIC определяется следующим образом:

\[AIC = -2 \log(\hat{L}) + 2 d\]

где \(\hat{L}\) — максимальное правдоподобие модели, а \(d\) — число параметров (также называемое числом степеней свободы в предыдущем разделе).

В определении BIC константа \(2\) заменяется на \(\log(N)\):

\[BIC = -2 \log(\hat{L}) + \log(N) d\]

где \(N\) — число образцов.

Для линейной гауссовской модели максимальное логарифмическое правдоподобие определяется следующим образом:

\[\log(\hat{L}) = - \frac{n}{2} \log(2 \pi) - \frac{n}{2} \log(\sigma^2) - \frac{\sum_{i=1}^{n} (y_i - \hat{y}_i)^2}{2\sigma^2}\]

где \(\sigma^2\) — оценка дисперсии шума, \(y_i\) и \(\hat{y}_i\) — соответственно истинные и предсказанные целевые значения, а \(n\) — число образцов.

Подстановка максимального логарифмического правдоподобия в формулу AIC даёт:

\[AIC = n \log(2 \pi \sigma^2) + \frac{\sum_{i=1}^{n} (y_i - \hat{y}_i)^2}{\sigma^2} + 2 d\]

Первый член этого выражения иногда отбрасывают, поскольку он является константой, если задано значение \(\sigma^2\). Кроме того, иногда утверждается, что AIC эквивалентен статистике \(C_p\) [12]. Однако в строгом смысле они эквивалентны лишь с точностью до константы и множителя.

Наконец, выше мы упоминали, что \(\sigma^2\) — это оценка дисперсии шума. В LassoLarsIC, если параметр noise_variance не задан (значение по умолчанию), дисперсия шума оценивается с помощью несмещённой оценки [13], определяемой следующим образом:

\[\sigma^2 = \frac{\sum_{i=1}^{n} (y_i - \hat{y}_i)^2}{n - p}\]

где \(p\) — число признаков, а \(\hat{y}_i\) — предсказанное целевое значение, полученное с помощью регрессии методом наименьших квадратов. Обратите внимание, что эта формула верна только при n_samples > n_features.

Литература

1.1.3.2.4. Сравнение с параметром регуляризации SVM

Соотношение между alpha и параметром регуляризации SVM, C, определяется как alpha = 1 / C или alpha = 1 / (n_samples * C) в зависимости от оценивателя и точной целевой функции, оптимизируемой моделью.

1.1.4. Многозадачное Лассо

MultiTaskLasso — это линейная модель, которая совместно оценивает разреженные коэффициенты для нескольких задач регрессии: y — это двумерный массив формы (n_samples, n_tasks). Ограничение заключается в том, что выбранные признаки одинаковы для всех задач регрессии, которые также называют задачами.

На следующем рисунке сравнивается расположение ненулевых элементов в матрице коэффициентов W, полученной с помощью обычного Лассо и MultiTaskLasso. Оценки Лассо дают разбросанные ненулевые элементы, тогда как ненулевые элементы MultiTaskLasso образуют целые столбцы.

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

Примеры

• Совместный отбор признаков с помощью многозадачного Лассо

1.1.5. Эластичная сеть

ElasticNet — это модель линейной регрессии, обучаемая с регуляризацией коэффициентов нормами \(\ell_1\) и \(\ell_2\). Такое сочетание позволяет обучать разреженную модель, в которой лишь немногие веса отличны от нуля, подобно Lasso, сохраняя при этом свойства регуляризации Ridge. Выпуклая комбинация норм \(\ell_1\) и \(\ell_2\) регулируется параметром l1_ratio.

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

Практическое преимущество компромисса между Лассо и Ridge заключается в том, что эластичная сеть наследует часть устойчивости Ridge к поворотам.

В этом случае минимизируемая целевая функция имеет вид

Класс ElasticNetCV можно использовать для настройки параметров alpha (\(\alpha\)) и l1_ratio (\(\rho\)) с помощью кросс-валидации.

Примеры

• Модели на основе L1 для разреженных сигналов
• Пути Лассо, Lasso-LARS и эластичной сети
• Подгонка эластичной сети с предварительно вычисленной матрицей Грама и взвешенными образцами

1.1.6. Многозадачная эластичная сеть

MultiTaskElasticNet — это модель эластичной сети, которая совместно оценивает разреженные коэффициенты для нескольких задач регрессии: Y — это двумерный массив формы (n_samples, n_tasks). Ограничение заключается в том, что выбранные признаки одинаковы для всех задач регрессии, которые также называют задачами.

Математически это линейная модель, обучаемая с регуляризацией смешанной нормой \(\ell_1\) \(\ell_2\) и нормой \(\ell_2\). Минимизируемая целевая функция имеет вид:

В реализации класса MultiTaskElasticNet для подгонки коэффициентов используется алгоритм покоординатного спуска.

Класс MultiTaskElasticNetCV можно использовать для настройки параметров alpha (\(\alpha\)) и l1_ratio (\(\rho\)) с помощью кросс-валидации.

1.1.7. Регрессия методом наименьшего угла

Регрессия методом наименьшего угла (LARS) — это алгоритм регрессии для данных высокой размерности, разработанный Брэдли Эфроном, Тревором Хасти, Иэном Джонстоном и Робертом Тибширани. LARS похожа на пошаговую регрессию с прямым отбором. На каждом шаге алгоритм находит признак, наиболее коррелирующий с целевой переменной. Если несколько признаков имеют одинаковую корреляцию, алгоритм вместо продолжения работы с тем же признаком движется в равноугольном направлении между признаками.

Преимущества LARS:

• Она эффективна с точки зрения вычислений в ситуациях, когда число признаков значительно превышает число выборок.
• Она работает так же быстро, как прямой отбор, и имеет тот же порядок сложности, что и обычный метод наименьших квадратов.
• Она строит полный кусочно-линейный путь решения, что полезно при кросс-валидации и других попытках настройки модели.
• Если два признака почти одинаково коррелируют с целевой переменной, их коэффициенты должны увеличиваться приблизительно с одинаковой скоростью. Таким образом, алгоритм ведет себя ожидаемым образом и является более устойчивым.
• Ее легко модифицировать для получения решений с помощью других оценивателей, например Lasso.

К недостаткам метода LARS относятся:

• Поскольку LARS основана на итеративном повторном подборе остатков, она, по-видимому, особенно чувствительна к влиянию шума. Эта проблема подробно рассматривается Вайсбергом в разделе обсуждения статьи Эфрона и др. (2004) в Annals of Statistics.

Модель LARS можно использовать с помощью оценивателя Lars или его низкоуровневой реализации lars_path или lars_path_gram.

1.1.8. LARS Lasso

LassoLars — это модель Lasso, реализованная с помощью алгоритма LARS. В отличие от реализации на основе покоординатного спуска, она дает точное решение, которое является кусочно-линейной функцией нормы коэффициентов.

>>> from sklearn import linear_model
>>> reg = linear_model.LassoLars(alpha=.1)
>>> reg.fit([[0, 0], [1, 1]], [0, 1])
LassoLars(alpha=0.1)
>>> reg.coef_
array([0.6, 0.        ])

Примеры

• Пути Lasso, Lasso-LARS и Elastic Net

Алгоритм LARS практически бесплатно строит полный путь коэффициентов в зависимости от параметра регуляризации, поэтому часто требуется получить этот путь с помощью одной из функций lars_path или lars_path_gram.

1.1.9. Ортогональный метод согласования (OMP)

OrthogonalMatchingPursuit и orthogonal_mp реализуют алгоритм OMP для аппроксимации подгонки линейной модели с ограничением на число ненулевых коэффициентов (то есть псевдонорму \(\ell_0\)).

Будучи методом прямого отбора признаков, подобно регрессии методом наименьшего угла, ортогональный метод согласования может аппроксимировать оптимальный вектор решения с фиксированным числом ненулевых элементов:

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

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

Примеры

• Ортогональный метод согласования

1.1.10. Байесовская регрессия

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

Это можно сделать, задав неинформативные априорные распределения для гиперпараметров модели. Регуляризация \(\ell_{2}\), используемая в гребневой регрессии и классификации, эквивалентна поиску оценки максимума апостериорной вероятности при гауссовском априорном распределении коэффициентов \(w\) с точностью \(\lambda^{-1}\). Вместо того чтобы вручную задавать lambda, его можно рассматривать как случайную величину, оцениваемую по данным.

Для получения полностью вероятностной модели предполагается, что выходная переменная \(y\) имеет гауссовское распределение с математическим ожиданием \(X w\):

где \(\alpha\) также рассматривается как случайная величина, оцениваемая по данным.

Преимущества байесовской регрессии:

• Она адаптируется к имеющимся данным.
• Ее можно использовать для включения параметров регуляризации в процедуру оценивания.

К недостаткам байесовской регрессии относятся:

• Вывод модели может занимать много времени.

>>> from sklearn import linear_model
>>> X = [[0., 0.], [1., 1.], [2., 2.], [3., 3.]]
>>> Y = [0., 1., 2., 3.]
>>> reg = linear_model.BayesianRidge()
>>> reg.fit(X, Y)
BayesianRidge()

>>> reg.predict([[1, 0.]])
array([0.50000013])

>>> reg.coef_
array([0.49999993, 0.49999993])

1.1.11. Логистическая регрессия

Логистическая регрессия реализована в LogisticRegression. Несмотря на своё название, в терминологии scikit-learn/машинного обучения она реализована как линейная модель классификации, а не регрессии. В литературе логистическая регрессия также известна как logit-регрессия, классификация с максимальной энтропией (MaxEnt) или лог-линейный классификатор. В этой модели вероятности возможных исходов отдельного испытания моделируются с помощью логистической функции.

Эта реализация позволяет обучать бинарную логистическую регрессию, логистическую регрессию по схеме «один против остальных» или мультиномиальную логистическую регрессию с необязательной регуляризацией \(\ell_1\), \(\ell_2\) или Elastic-Net.

Примеры

• Штраф L1 и разреженность в логистической регрессии
• Траектория регуляризации логистической регрессии с L1
• Границы принятия решений для мультиномиальной логистической регрессии и регрессии по схеме «один против остальных»
• Многоклассовая разреженная логистическая регрессия для 20newgroups
• Классификация MNIST с использованием мультиномиальной логистической регрессии и L1
• Визуализация вероятности классификации

1.1.12. Обобщённые линейные модели

Обобщённые линейные модели (GLM) расширяют линейные модели двумя способами [10]. Во-первых, прогнозируемые значения \(\hat{y}\) связываются с линейной комбинацией входных переменных \(X\) с помощью обратной функции связи \(h\) следующим образом

Во-вторых, функция квадратичной ошибки заменяется единичной девиансой \(d\) распределения из экспоненциального семейства (или, точнее, воспроизводящей экспоненциальной дисперсионной модели (EDM) [11]).

Задача минимизации принимает вид:

где \(\alpha\) — штраф за L2-регуляризацию. При наличии весов выборок среднее становится взвешенным средним.

В следующей таблице приведены некоторые конкретные EDM и их единичные девиансы:

Функции плотности вероятности (PDF) этих распределений показаны на следующем рисунке,

Распределение Бернулли — это дискретное распределение вероятностей, моделирующее испытание Бернулли — событие, у которого есть только два взаимоисключающих исхода. Категориальное распределение является обобщением распределения Бернулли для категориальной случайной величины. Если случайная величина в распределении Бернулли имеет два возможных исхода, то категориальная случайная величина может принимать одну из K возможных категорий, причём вероятность каждой категории задаётся отдельно.

Выбор распределения зависит от решаемой задачи:

• Если целевые значения \(y\) — это количества (неотрицательные целочисленные значения) или относительные частоты (неотрицательные значения), можно использовать распределение Пуассона с логарифмической функцией связи.
• Если целевые значения положительны и имеют асимметричное распределение, можно попробовать гамма-распределение с логарифмической функцией связи.
• Если целевые значения, по-видимому, имеют более тяжёлые хвосты, чем гамма-распределение, можно попробовать обратное гауссово распределение (или даже распределения семейства Твиди с более высокими степенями дисперсии).
• Если целевые значения \(y\) — это вероятности, можно использовать распределение Бернулли. Распределение Бернулли с логит-функцией связи можно использовать для бинарной классификации. Категориальное распределение с функцией связи softmax можно использовать для многоклассовой классификации.

Ссылки

>>> from sklearn.linear_model import TweedieRegressor
>>> reg = TweedieRegressor(power=1, alpha=0.5, link='log')
>>> reg.fit([[0, 0], [0, 1], [2, 2]], [0, 1, 2])
TweedieRegressor(alpha=0.5, link='log', power=1)
>>> reg.coef_
array([0.2463, 0.4337])
>>> reg.intercept_
np.float64(-0.7638)

1.1.13. Стохастический градиентный спуск — SGD

Стохастический градиентный спуск — простой, но очень эффективный способ обучения линейных моделей. Он особенно полезен, когда число выборок (и число признаков) очень велико. Метод partial_fit позволяет обучаться в онлайн-режиме и обрабатывать данные, не помещающиеся в память.

Классы SGDClassifier и SGDRegressor позволяют обучать линейные модели для классификации и регрессии с использованием различных выпуклых функций потерь и различных штрафов. Например, с параметром loss="log" SGDClassifier обучает модель логистической регрессии, а с параметром loss="hinge" — линейную машину опорных векторов (SVM).

Дополнительные сведения см. в специальном разделе документации «Стохастический градиентный спуск».

1.1.14. Робастная регрессия: выбросы и ошибки моделирования

Робастная регрессия предназначена для построения регрессионной модели при наличии искажённых данных: выбросов или ошибок в модели.

1.1.14.2. RANSAC: консенсус случайной выборки

RANSAC (RANdom SAmple Consensus) строит модель на основе случайных подмножеств точек без выбросов из полного набора данных.

RANSAC — недетерминированный алгоритм, который с некоторой вероятностью выдаёт приемлемый результат; эта вероятность зависит от числа итераций (см. параметр max_trials). Обычно его используют в задачах линейной и нелинейной регрессии; он особенно популярен в фотограмметрическом компьютерном зрении.

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

Примеры

• Оценивание робастной линейной модели с помощью RANSAC
• Подгонка робастного линейного оценивателя

Описание алгоритма

На каждой итерации выполняются следующие шаги:

1. Выбрать min_samples случайных образцов из исходных данных и проверить допустимость набора данных (см. is_data_valid).
2. Подогнать модель по случайному подмножеству (estimator.fit) и проверить допустимость оценённой модели (см. is_model_valid).
3. Классифицировать все данные как точки без выбросов или выбросы, вычислив остатки относительно оценённой модели (estimator.predict(X) - y): все образцы данных с абсолютными остатками, меньшими или равными residual_threshold, считаются точками без выбросов.
4. Сохранить подогнанную модель как наилучшую, если число образцов без выбросов максимально. Если в текущей оценённой модели столько же точек без выбросов, она считается наилучшей только при более высокой оценке.

Эти шаги выполняются либо максимальное число раз (max_trials), либо до выполнения одного из специальных критериев остановки (см. stop_n_inliers и stop_score). Итоговая модель оценивается по всем образцам без выбросов (консенсусному набору) наилучшей из ранее определённых моделей.

Функции is_data_valid и is_model_valid позволяют выявлять и отбрасывать вырожденные комбинации случайных подвыборок. Если оценённая модель не нужна для выявления вырожденных случаев, следует использовать is_data_valid: она вызывается до подгонки модели и тем самым повышает вычислительную эффективность.

Ссылки

• https://en.wikipedia.org/wiki/RANSAC
• «Консенсус случайной выборки: парадигма подгонки моделей с применением к анализу изображений и автоматизированной картографии» Martin A. Fischler и Robert C. Bolles — SRI International (1981)
• «Оценка эффективности семейства RANSAC» Sunglok Choi, Taemin Kim и Wonpil Yu — BMVC (2009)

1.1.14.3. Оцениватель Theil–Sen: оцениватель на основе обобщённой медианы

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

Примеры

• Регрессия Theil–Sen
• Подгонка робастного линейного оценивателя

Теоретические замечания

По асимптотической эффективности и несмещённости оцениватель TheilSenRegressor сопоставим с методом обыкновенных наименьших квадратов (OLS). В отличие от OLS, Theil–Sen — непараметрический метод, то есть он не делает предположений о распределении данных. Поскольку Theil–Sen — оцениватель на основе медианы, он более устойчив к искажённым данным, то есть выбросам. В одномерном случае точка разрушения Theil–Sen составляет около 29,3 % для простой линейной регрессии, а значит, он может выдерживать произвольное искажение до 29,3 % данных.

Реализация TheilSenRegressor в scikit-learn следует обобщению для многомерной линейной регрессионной модели [14], использующему пространственную медиану — обобщение медианы для нескольких измерений [15].

Вычислительная сложность Theil–Sen по времени и памяти растёт пропорционально

\[\binom{n_{\text{samples}}}{n_{\text{subsamples}}}\]

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

Ссылки

См. также страницу Википедии

1.1.14.4. Регрессия Хьюбера

HuberRegressor отличается от Ridge тем, что применяет линейную функцию потерь к образцам, определённым как выбросы параметром epsilon. Образец классифицируется как точка без выброса, если абсолютная ошибка для него меньше порога epsilon. От оценивателей TheilSenRegressor и RANSACRegressor он отличается тем, что не игнорирует влияние выбросов, а лишь присваивает им меньший вес.

Примеры

• HuberRegressor и Ridge на наборе данных с сильными выбросами

Математические подробности

HuberRegressor минимизирует

\[\min_{w, \sigma} {\sum_{i=1}^n\left(\sigma + H_{\epsilon}\left(\frac{X_{i}w - y_{i}}{\sigma}\right)\sigma\right) + \alpha {||w||_2}^2}\]

где функция потерь задаётся выражением

\[\begin{split}H_{\epsilon}(z) = \begin{cases} z^2, & \text {if } |z| < \epsilon, \\ 2\epsilon|z| - \epsilon^2, & \text{otherwise} \end{cases}\end{split}\]

Для достижения статистической эффективности 95 % рекомендуется задать параметру epsilon значение 1.35.

Ссылки

• Peter J. Huber, Elvezio M. Ronchetti: «Робастная статистика», оценки сопутствующего масштаба, с. 172.

HuberRegressor отличается от использования SGDRegressor с параметром loss, равным huber, следующим образом.

• HuberRegressor инвариантен к масштабированию. Если задать epsilon, масштабирование X и y в большую или меньшую сторону с разными коэффициентами сохранит прежнюю устойчивость к выбросам. В отличие от него, в SGDRegressor параметр epsilon необходимо задавать заново при масштабировании X и y.
• HuberRegressor должен быть эффективнее на данных с небольшим числом образцов, тогда как SGDRegressor требуется несколько проходов по обучающим данным, чтобы обеспечить такую же робастность.

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

1.1.15. Квантильная регрессия

Квантильная регрессия оценивает медиану или другие квантили \(y\) при условии \(X\), тогда как обычный метод наименьших квадратов (OLS) оценивает условное среднее.

Квантильная регрессия может быть полезна, если требуется предсказать интервал, а не отдельное значение. Иногда интервалы предсказания вычисляются исходя из предположения, что ошибка предсказания распределена нормально с нулевым средним и постоянной дисперсией. Квантильная регрессия даёт разумные интервалы предсказания даже при непостоянной (но предсказуемой) дисперсии ошибок или ненормальном распределении.

Условные квантили можно оценивать и моделями, отличными от линейных, минимизируя pinball loss. Например, GradientBoostingRegressor может предсказывать условные квантили, если параметру loss присвоено значение "quantile", а параметру alpha — значение квантили, которую требуется предсказать. См. пример Интервалы предсказания для регрессии с градиентным бустингом.

Большинство реализаций квантильной регрессии основано на задаче линейного программирования. Текущая реализация использует scipy.optimize.linprog.

Примеры

• Квантильная регрессия

1.1.16. Полиномиальная регрессия: расширение линейных моделей с помощью базисных функций

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

Вот пример применения этой идеи к одномерным данным с использованием полиномиальных признаков различной степени:

Этот рисунок создан с помощью преобразователя PolynomialFeatures, который преобразует входную матрицу данных в новую матрицу заданной степени. Его можно использовать следующим образом:

>>> from sklearn.preprocessing import PolynomialFeatures
>>> import numpy as np
>>> X = np.arange(6).reshape(3, 2)
>>> X
array([[0, 1],
       [2, 3],
       [4, 5]])
>>> poly = PolynomialFeatures(degree=2)
>>> poly.fit_transform(X)
array([[ 1.,  0.,  1.,  0.,  0.,  1.],
       [ 1.,  2.,  3.,  4.,  6.,  9.],
       [ 1.,  4.,  5., 16., 20., 25.]])

Признаки X преобразованы из \([x_1, x_2]\) в \([1, x_1, x_2, x_1^2, x_1 x_2, x_2^2]\), и теперь их можно использовать в любой линейной модели.

Такую предварительную обработку можно упростить с помощью инструментов Pipeline. Один объект, представляющий простую полиномиальную регрессию, можно создать и использовать следующим образом:

>>> from sklearn.preprocessing import PolynomialFeatures
>>> from sklearn.linear_model import LinearRegression
>>> from sklearn.pipeline import Pipeline
>>> import numpy as np
>>> model = Pipeline([('poly', PolynomialFeatures(degree=3)),
...                   ('linear', LinearRegression(fit_intercept=False))])
>>> # fit to an order-3 polynomial data
>>> x = np.arange(5)
>>> y = 3 - 2 * x + x ** 2 - x ** 3
>>> model = model.fit(x[:, np.newaxis], y)
>>> model.named_steps['linear'].coef_
array([ 3., -2.,  1., -1.])

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

В некоторых случаях необязательно включать высшие степени отдельных признаков; достаточно так называемых признаков взаимодействия, которые перемножают не более \(d\) различных признаков. Их можно получить с помощью PolynomialFeatures, задав параметр interaction_only=True.

Например, для булевых признаков \(x_i^n = x_i\) при любом \(n\), поэтому такие степени бесполезны; однако \(x_i x_j\) представляет собой конъюнкцию двух булевых значений. Таким образом, мы можем решить задачу XOR с помощью линейного классификатора:

>>> from sklearn.linear_model import Perceptron
>>> from sklearn.preprocessing import PolynomialFeatures
>>> import numpy as np
>>> X = np.array([[0, 0], [0, 1], [1, 0], [1, 1]])
>>> y = X[:, 0] ^ X[:, 1]
>>> y
array([0, 1, 1, 0])
>>> X = PolynomialFeatures(interaction_only=True).fit_transform(X).astype(int)
>>> X
array([[1, 0, 0, 0],
       [1, 0, 1, 0],
       [1, 1, 0, 0],
       [1, 1, 1, 1]])
>>> clf = Perceptron(fit_intercept=False, max_iter=10, tol=None,
...                  shuffle=False).fit(X, y)

И классификатор выдаёт идеальные «предсказания»:

>>> clf.predict(X)
array([0, 1, 1, 0])
>>> clf.score(X, y)
1.0