Меню

Линейная регрессия среднеквадратичная ошибка

Среднеквадратичная ошибка (Mean Squared Error) – Среднее арифметическое (Mean) квадратов разностей между предсказанными и реальными значениями Модели (Model) Машинного обучения (ML):

MSE как среднее дистанций между предсказаниями и реальными наблюдениями

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

$$MSE = frac{1}{n} × sum_{i=1}^n (y_i — widetilde{y}_i)^2$$
$$MSEspace{}{–}space{Среднеквадратическая}space{ошибка,}$$
$$nspace{}{–}space{количество}space{наблюдений,}$$
$$y_ispace{}{–}space{фактическая}space{координата}space{наблюдения,}$$
$$widetilde{y}_ispace{}{–}space{предсказанная}space{координата}space{наблюдения,}$$

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

Пример. Исследуем линейную регрессию, изображенную на графике выше, и установим величину среднеквадратической Ошибки (Error). Фактические координаты точек-Наблюдений (Observation) выглядят следующим образом:

Мы имеем дело с Линейной регрессией (Linear Regression), потому уравнение, предсказывающее положение записей, можно представить с помощью формулы:

$$y = M * x + b$$
$$yspace{–}space{значение}space{координаты}space{оси}space{y,}$$
$$Mspace{–}space{уклон}space{прямой}$$
$$xspace{–}space{значение}space{координаты}space{оси}space{x,}$$
$$bspace{–}space{смещение}space{прямой}space{относительно}space{начала}space{координат}$$

Параметры M и b уравнения нам, к счастью, известны в данном обучающем примере, и потому уравнение выглядит следующим образом:

$$y = 0,5252 * x + 17,306$$

Зная координаты реальных записей и уравнение линейной регрессии, мы можем восстановить полные координаты предсказанных наблюдений, обозначенных серыми точками на графике выше. Простой подстановкой значения координаты x в уравнение мы рассчитаем значение координаты ỹ:

Рассчитаем квадрат разницы между Y и Ỹ:

Сумма таких квадратов равна 4 445. Осталось только разделить это число на количество наблюдений (9):

$$MSE = frac{1}{9} × 4445 = 493$$

Само по себе число в такой ситуации становится показательным, когда Дата-сайентист (Data Scientist) предпринимает попытки улучшить предсказательную способность модели и сравнивает MSE каждой итерации, выбирая такое уравнение, что сгенерирует наименьшую погрешность в предсказаниях.

MSE и Scikit-learn

Среднеквадратическую ошибку можно вычислить с помощью SkLearn. Для начала импортируем функцию:

import sklearn
from sklearn.metrics import mean_squared_error

Инициализируем крошечные списки, содержащие реальные и предсказанные координаты y:

y_true = [5, 41, 70, 77, 134, 68, 138, 101, 131]
y_pred = [23, 35, 55, 90, 93, 103, 118, 121, 129]

Инициируем функцию mean_squared_error(), которая рассчитает MSE тем же способом, что и формула выше:

mean_squared_error(y_true, y_pred)

Интересно, что конечный результат на 3 отличается от расчетов с помощью Apple Numbers:

496.0

Ноутбук, не требующий дополнительной настройки на момент написания статьи, можно скачать здесь.

Автор оригинальной статьи: @mmoshikoo

Фото: @tobyelliott

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

«Хорошая» аналитическая модель должна удовлетворять двум, зачастую противоречивым, требованиям — как можно лучше соответствовать данным и при этом быть удобной для интерпретации пользователем. Действительно, повышение соответствия модели данным как правило связано с её усложнением (в случае регрессии — увеличением числа входных переменных модели). А чем сложнее модель, тем ниже её интерпретируемость.

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

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

  • Среднеквадратичная ошибка (MSE).
  • Корень из среднеквадратичной ошибки (RMSE).
  • Среднеквадратичная ошибка в процентах (MSPE).
  • Средняя абсолютная ошибка (MAE).
  • Средняя абсолютная ошибка в процентах (MAPE).
  • Cимметричная средняя абсолютная процентная ошибка (SMAPE).
  • Средняя абсолютная масштабированная ошибка (MASE)
  • Средняя относительная ошибка (MRE).
  • Среднеквадратичная логарифмическая ошибка (RMSLE).
  • Коэффициент детерминации R-квадрат.
  • Скорректированный коэффициент детеминации.

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

Рисунок 1. Линейная регрессия

Наклонная прямая представляет собой линию регрессии с переменной, на которой расположены точки, соответствующие предсказанным значениям выходной переменной widehat{y} (кружки синего цвета). Оранжевые кружки представляют фактические (наблюдаемые) значения y . Расстояния между ними и линией регрессии — это ошибка предсказания модели y-widehat{y} (невязка, остатки). Именно с её использованием вычисляются все приведённые в статье меры качества.

Горизонтальная линия представляет собой модель простого среднего, где коэффициент при независимой переменной x равен нулю, и остаётся только свободный член b, который становится равным среднему арифметическому фактических значений выходной переменной, т.е. b=overline{y}. Очевидно, что такая модель для любого значения входной переменной будет выдавать одно и то же значение выходной — overline{y}.

В линейной регрессии такая модель рассматривается как «бесполезная», хуже которой работает только «случайный угадыватель». Однако, она используется для оценки, насколько дисперсия фактических значений y относительно линии среднего, больше, чем относительно линии регрессии с переменной, т.е. насколько модель с переменной лучше «бесполезной».

MSE

Среднеквадратичная ошибка (Mean Squared Error) применяется в случаях, когда требуется подчеркнуть большие ошибки и выбрать модель, которая дает меньше именно больших ошибок. Большие значения ошибок становятся заметнее за счет квадратичной зависимости.

Действительно, допустим модель допустила на двух примерах ошибки 5 и 10. В абсолютном выражении они отличаются в два раза, но если их возвести в квадрат, получив 25 и 100 соответственно, то отличие будет уже в четыре раза. Таким образом модель, которая обеспечивает меньшее значение MSE допускает меньше именно больших ошибок.

MSE рассчитывается по формуле:

MSE=frac{1}{n}sumlimits_{i=1}^{n}(y_{i}-widehat{y}_{i})^{2},

где n — количество наблюдений по которым строится модель и количество прогнозов, y_{i} — фактические значение зависимой переменной для i-го наблюдения, widehat{y}_{i} — значение зависимой переменной, предсказанное моделью.

Таким образом, можно сделать вывод, что MSE настроена на отражение влияния именно больших ошибок на качество модели.

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

RMSE

Корень из среднеквадратичной ошибки (Root Mean Squared Error) вычисляется просто как квадратный корень из MSE:

RMSE=sqrt{frac{1}{n}sumlimits_{i=1}^{n}(y_{i}-widehat{y_{i}})^{2}}

MSE и RMSE могут минимизироваться с помощью одного и того же функционала, поскольку квадратный корень является неубывающей функцией. Например, если у нас есть два набора результатов работы модели, A и B, и MSE для A больше, чем MSE для B, то мы можем быть уверены, что RMSE для A больше RMSE для B. Справедливо и обратное: если MSE(A)<MSE(B), то и RMSE(A)<RMSE(B).

Следовательно, сравнение моделей с помощью RMSE даст такой же результат, что и для MSE. Однако с MSE работать несколько проще, поэтому она более популярна у аналитиков. Кроме этого, имеется небольшая разница между этими двумя ошибками при оптимизации с использованием градиента:

frac{partial RMSE}{partial widehat{y}_{i}}=frac{1}{2sqrt{MSE}}frac{partial MSE}{partial widehat{y}_{i}}

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

Влияние каждой ошибки на RMSE пропорционально величине квадрата ошибки. Поэтому большие ошибки оказывают непропорционально большое влияние на RMSE. Следовательно, RMSE можно считать чувствительной к аномальным значениям.

MSPE

Среднеквадратичная ошибка в процентах (Mean Squared Percentage Error) представляет собой относительную ошибку, где разность между наблюдаемым и фактическим значениями делится на наблюдаемое значение и выражается в процентах:

MSPE=frac{100}{n}sumlimits_{i=1}^{n}left ( frac{y_{i}-widehat{y}_{i}}{y_{i}} right )^{2}

Проблемой при использовании MSPE является то, что, если наблюдаемое значение выходной переменной равно 0, значение ошибки становится неопределённым.

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

MAE

Cредняя абсолютная ошибка (Mean Absolute Error) вычисляется следующим образом:

MAE=frac{1}{n}sumlimits_{i=1}^{n}left | y_{i}-widehat{y}_{i} right |

Т.е. MAE рассчитывается как среднее абсолютных разностей между наблюдаемым и предсказанным значениями. В отличие от MSE и RMSE она является линейной оценкой, а это значит, что все ошибки в среднем взвешены одинаково. Например, разница между 0 и 10 будет вдвое больше разницы между 0 и 5. Для MSE и RMSE, как отмечено выше, это не так.

Поэтому MAE широко используется, например, в финансовой сфере, где ошибка в 10 долларов должна интерпретироваться как в два раза худшая, чем ошибка в 5 долларов.

MAPE

Средняя абсолютная процентная ошибка (Mean Absolute Percentage Error) вычисляется следующим образом:

MAPE=frac{100}{n}sumlimits_{i=1}^{n}frac{left | y_{i}-widehat{y_{i}} right |}{left | y_{i} right |}

Эта ошибка не имеет размерности и очень проста в интерпретации. Её можно выражать как в долях, так и в процентах. Если получилось, например, что MAPE=11.4, то это говорит о том, что ошибка составила 11.4% от фактического значения.

SMAPE

Cимметричная средняя абсолютная процентная ошибка (Symmetric Mean Absolute Percentage Error) — это мера точности, основанная на процентных (или относительных) ошибках. Обычно определяется следующим образом:

SMAPE=frac{100}{n}sumlimits_{i=1}^{n}frac{left | y_{i}-widehat{y_{i}} right |}{(left | y_{i} right |+left | widehat{y}_{i} right |)/2}

Т.е. абсолютная разность между наблюдаемым и предсказанным значениями делится на полусумму их модулей. В отличие от обычной MAPE, симметричная имеет ограничение на диапазон значений. В приведённой формуле он составляет от 0 до 200%. Однако, поскольку диапазон от 0 до 100% гораздо удобнее интерпретировать, часто используют формулу, где отсутствует деление знаменателя на 2.

Одной из возможных проблем SMAPE является неполная симметрия, поскольку в разных диапазонах ошибка вычисляется неодинаково. Это иллюстрируется следующим примером: если y_{i}=100 и widehat{y}_{i}=110, то SMAPE=4.76, а если y_{i}=100 и widehat{y}_{i}=90, то SMAPE=5.26.

Ограничение SMAPE заключается в том, что, если наблюдаемое или предсказанное значение равно 0, ошибка резко возрастет до верхнего предела (200% или 100%).

MASE

Средняя абсолютная масштабированная ошибка (Mean absolute scaled error) — это показатель, который позволяет сравнивать две модели. Если поместить MAE для новой модели в числитель, а MAE для исходной модели в знаменатель, то полученное отношение и будет равно MASE. Если значение MASE меньше 1, то новая модель работает лучше, если MASE равно 1, то модели работают одинаково, а если значение MASE больше 1, то исходная модель работает лучше, чем новая модель. Формула для расчета MASE имеет вид:

MASE=frac{MAE_{i}}{MAE_{j}}

MASE симметрична и устойчива к выбросам.

MRE

Средняя относительная ошибка (Mean Relative Error) вычисляется по формуле:

MRE=frac{1}{n}sumlimits_{i=1}^{n}frac{left | y_{i}-widehat{y}_{i}right |}{left | y_{i} right |}

Несложно увидеть, что данная мера показывает величину абсолютной ошибки относительно фактического значения выходной переменной (поэтому иногда эту ошибку называют также средней относительной абсолютной ошибкой, MRAE). Действительно, если значение абсолютной ошибки, скажем, равно 10, то сложно сказать много это или мало. Например, относительно значения выходной переменной, равного 20, это составляет 50%, что достаточно много. Однако относительно значения выходной переменной, равного 100, это будет уже 10%, что является вполне нормальным результатом.

Очевидно, что при вычислении MRE нельзя применять наблюдения, в которых y_{i}=0.

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

RMSLE

Среднеквадратичная логарифмическая ошибка (Root Mean Squared Logarithmic Error) представляет собой RMSE, вычисленную в логарифмическом масштабе:

RMSLE=sqrt{frac{1}{n}sumlimits_{i=1}^{n}(log(widehat{y}_{i}+1)-log{(y_{i}+1}))^{2}}

Константы, равные 1, добавляемые в скобках, необходимы чтобы не допустить обращения в 0 выражения под логарифмом, поскольку логарифм нуля не существует.

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

R-квадрат

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

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

Именно такой мерой и является коэффициент детерминации (Coefficient of determination), который показывает долю дисперсии зависимой переменной, объяснённой с помощью регрессионной модели. Наиболее общей формулой для вычисления коэффициента детерминации является следующая:

R^{2}=1-frac{sumlimits_{i=1}^{n}(widehat{y}_{i}-y_{i})^{2}}{sumlimits_{i=1}^{n}({overline{y}}_{i}-y_{i})^{2}}

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

Главным преимуществом коэффициента детерминации перед мерами, основанными на ошибках, является его инвариантность к масштабу данных. Кроме того, он всегда изменяется в диапазоне от −∞ до 1. При этом значения близкие к 1 указывают на высокую степень соответствия модели данным. Очевидно, что это имеет место, когда отношение в формуле стремится к 0, т.е. ошибка модели с переменными намного меньше ошибки модели с константой. R^{2}=0 показывает, что между независимой и зависимой переменными модели имеет место функциональная зависимость.

Когда значение коэффициента близко к 0 (т.е. ошибка модели с переменными примерно равна ошибке модели только с константой), это указывает на низкое соответствие модели данным, когда модель с переменными работает не лучше модели с константой.

Кроме этого, бывают ситуации, когда коэффициент R^{2} принимает отрицательные значения (обычно небольшие). Это произойдёт, если ошибка модели среднего становится меньше ошибки модели с переменной. В этом случае оказывается, что добавление в модель с константой некоторой переменной только ухудшает её (т.е. регрессионная модель с переменной работает хуже, чем предсказание с помощью простой средней).

На практике используют следующую шкалу оценок. Модель, для которой R^{2}>0.5, является удовлетворительной. Если R^{2}>0.8, то модель рассматривается как очень хорошая. Значения, меньшие 0.5 говорят о том, что модель плохая.

Скорректированный R-квадрат

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

В связи с этим возникают две проблемы. Первая заключается в том, что не все переменные, добавляемые в модель, могут значимо увеличивать её точность, но при этом всегда увеличивают её сложность. Вторая проблема — с помощью коэффициента детерминации нельзя сравнивать модели с разным числом переменных. Чтобы преодолеть эти проблемы используют альтернативные показатели, одним из которых является скорректированный коэффициент детерминации (Adjasted coefficient of determinftion).

Скорректированный коэффициент детерминации даёт возможность сравнивать модели с разным числом переменных так, чтобы их число не влияло на статистику R^{2}, и накладывает штраф за дополнительно включённые в модель переменные. Вычисляется по формуле:

R_{adj}^{2}=1-frac{sumlimits_{i=1}^{n}(widehat{y}_{i}-y_{i})^{2}/(n-k)}{sumlimits_{i=1}^{n}({overline{y}}_{i}-y_{i})^{2}/(n-1)}

где n — число наблюдений, на основе которых строится модель, k — количество переменных в модели.

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

Сравнение метрик

Резюмируем преимущества и недостатки каждой приведённой метрики в следующей таблице:

Мера Сильные стороны Слабые стороны
MSE Позволяет подчеркнуть большие отклонения, простота вычисления. Имеет тенденцию занижать качество модели, чувствительна к выбросам. Сложность интерпретации из-за квадратичной зависимости.
RMSE Простота интерпретации, поскольку измеряется в тех же единицах, что и целевая переменная. Имеет тенденцию занижать качество модели, чувствительна к выбросам.
MSPE Нечувствительна к выбросам. Хорошо интерпретируема, поскольку имеет линейный характер. Поскольку вклад всех ошибок отдельных наблюдений взвешивается одинаково, не позволяет подчёркивать большие и малые ошибки.
MAPE Является безразмерной величиной, поэтому её интерпретация не зависит от предметной области. Нельзя использовать для наблюдений, в которых значения выходной переменной равны нулю.
SMAPE Позволяет корректно работать с предсказанными значениями независимо от того больше они фактического, или меньше. Приближение к нулю фактического или предсказанного значения приводит к резкому росту ошибки, поскольку в знаменателе присутствует как фактическое, так и предсказанное значения.
MASE Не зависит от масштаба данных, является симметричной: положительные и отрицательные отклонения от фактического значения учитываются одинаково. Устойчива к выбросам. Позволяет сравнивать модели. Сложность интерпретации.
MRE Позволяет оценить величину ошибки относительно значения целевой переменной. Неприменима для наблюдений с нулевым значением выходной переменной.
RMSLE Логарифмирование позволяет сделать величину ошибки более устойчивой, когда разность между фактическим и предсказанным значениями различается на порядок и выше Может быть затруднена интерпретация из-за нелинейности.
R-квадрат Универсальность, простота интерпретации. Возрастает даже при включении в модель бесполезных переменных. Плохо работает когда входные переменные зависимы.
R-квадрат скорр. Корректно отражает вклад каждой переменной в модель. Плохо работает, когда входные переменные зависимы.

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

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

Другие материалы по теме:

Отбор переменных в моделях линейной регрессии

Репрезентативность выборочных данных

Логистическая регрессия и ROC-анализ — математический аппарат

В предыдущих заметках предметом анализа часто становилась отдельная числовая переменная, например, доходность взаимных фондов, время загрузки Web-страницы или объем потребления безалкогольных напитков. В настоящей и следующих заметках мы рассмотрим методы предсказания значений числовой переменной в зависимости от значений одной или нескольких других числовых переменных. [1]

Материал будет проиллюстрирован сквозным примером. Прогнозирование объема продаж в магазине одежды. Сеть магазинов уцененной одежды Sunflowers на протяжении 25 лет постоянно расширялась. Однако в настоящее время у компании нет систематического подхода к выбору новых торговых точек. Место, в котором компания собирается открыть новый магазин, определяется на основе субъективных соображений. Критериями выбора являются выгодные условия аренды или представления менеджера об идеальном местоположении магазина. Представьте, что вы — руководитель отдела специальных проектов и планирования. Вам поручили разработать стратегический план открытия новых магазинов. Этот план должен содержать прогноз годового объема продаж во вновь открываемых магазинах. Вы полагаете, что торговая площадь непосредственно связана с объемом выручки, и хотите учесть этот факт в процессе принятия решения. Как разработать статистическую модель, позволяющую прогнозировать годовой объем продаж на основе размера нового магазина?

Как правило, для предсказания значений переменной используется регрессионный анализ. Его цель — разработать статистическую модель, позволяющую предсказывать значения зависимой переменной, или отклика, по значениям, по крайней мере одной, независимой, или объясняющей, переменной. В настоящей заметке мы рассмотрим простую линейную регрессию — статистический метод, позволяющий предсказывать значения зависимой переменной Y по значениям независимой переменной X. В последующих заметках будет описана модель множественной регрессии, предназначенная для предсказания значений независимой переменной Y по значениям нескольких зависимых переменных (Х1, Х2, …, Xk). [2]

Скачать заметку в формате Word или pdf, примеры в формате Excel2013

Виды регрессионных моделей

В заметке Представление числовых данных в виде таблиц и диаграмм для иллюстрации зависимости между переменными X и Y использовалась диаграмма разброса. На ней значения переменной X откладывались по горизонтальной оси, а значения переменной Y — по вертикальной. Зависимость между двумя переменными может быть разной: от самой простой до крайне сложной. Пример простейшей (линейной) зависимости показан на рис. 1.

Рис. 1. Положительная линейная зависимость

Простая линейная регрессия:

(1) Yi = β0 + β1Xi + εi

где β0 — сдвиг (длина отрезка, отсекаемого на координатной оси прямой Y), β1 — наклон прямой Y, εi— случайная ошибка переменной Y в i-м наблюдении.

В этой модели наклон β1 представляет собой количество единиц измерения переменной Y, приходящихся на одну единицу измерения переменной X. Эта величина характеризует среднюю величину изменения переменной Y (положительного или отрицательного) на заданном отрезке оси X. Сдвиг β0 представляет собой среднее значение переменной Y, когда переменная X равна 0. Последний компонент модели εi является случайной ошибкой переменной Y в i-м наблюдении. Выбор подходящей математической модели зависит от распределения значений переменных X и Y на диаграмме разброса. Различные виды зависимости переменных показаны на рис. 2.

Рис. 2. Диаграммы разброса, иллюстрирующие разные виды зависимостей

На панели А значения переменной Y почти линейно возрастают с увеличением переменной X. Этот рисунок аналогичен рис. 1, иллюстрирующему положительную зависимость между размером магазина (в квадратных футах) и годовым объемом продаж. Панель Б является примером отрицательной линейной зависимости. Если переменная X возрастает, переменная Y в целом убывает. Примером этой зависимости является связь между стоимостью конкретного товара и объемом продаж. На панели В показан набор данных, в котором переменные X и Y практически не зависят друг от друга. Каждому значению переменной X соответствуют как большие, так и малые значения переменной Y. Данные, приведенные на панели Г, демонстрируют криволинейную зависимость между переменными X и Y. Значения переменной Y возрастают при увеличении переменной X, однако скорость роста после определенных значений переменной X падает. Примером положительной криволинейной зависимости является связь между возрастом и стоимостью обслуживания автомобилей. По мере старения машины стоимость ее обслуживания сначала резко возрастает, однако после определенного уровня стабилизируется. Панель Д демонстрирует параболическую U-образную форму зависимости между переменными X и Y. По мере увеличения значений переменной X значения переменной Y сначала убывают, а затем возрастают. Примером такой зависимости является связь между количеством ошибок, совершенных за час работы, и количеством отработанных часов. Сначала работник осваивается и делает много ошибок, потом привыкает, и количество ошибок уменьшается, однако после определенного момента он начинает чувствовать усталость, и число ошибок увеличивается. На панели Е показана экспоненциальная зависимость между переменными X и Y. В этом случае переменная Y сначала очень быстро убывает при возрастании переменной X, однако скорость этого убывания постепенно падает. Например, стоимость автомобиля при перепродаже экспоненциально зависит от его возраста. Если перепродавать автомобиль в течение первого года, его цена резко падает, однако впоследствии ее падение постепенно замедляется.

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

Вывод уравнения простой линейной регрессии

Вернемся к сценарию, изложенному в начале главы. Наша цель — предсказать объем годовых продаж для всех новых магазинов, зная их размеры. Для оценки зависимости между размером магазина (в квадратных футах) и объемом его годовых продаж создадим выборки из 14 магазинов (рис. 3).

Рис. 3. Площади и годовые объемы продаж 14 магазинов сети Sunflowers: (а) исходные данные; (б) диаграмма разброса

Анализ рис. 3 показывает, что между площадью магазина X и годовым объемом продаж Y существует положительная зависимость. Если площадь магазина увеличивается, объем продаж возрастает почти линейно. Таким образом, наиболее подходящей для исследования является линейная модель. Остается лишь определить, какая из линейных моделей точнее остальных описывает зависимость между анализируемыми переменными.

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

Данные, представленные на рис. 1а, получены для случайной выборки магазинов. Если верны некоторые предположения (об этом чуть позже), в качестве оценки параметров генеральной совокупности (β0 и β1) можно использовать сдвиг b0 и наклон b1 прямой Y. Таким образом, уравнение простой линейной регрессии принимает следующий вид:

где — предсказанное значение переменной Y для i-гo наблюдения, Xi — значение переменной X в i-м наблюдении.

Для того чтобы предсказать значение переменной Y, в уравнении (2) необходимо определить два коэффициента регрессии — сдвиг b0 и наклон b1 прямой Y. Вычислив эти параметры, проведем прямую на диаграмме разброса. Затем исследователь может визуально оценить, насколько близка регрессионная прямая к точкам наблюдения. Простая линейная регрессия позволяет найти прямую линию, максимально приближенную к точкам наблюдения. Критерии соответствия можно задать разными способами. Возможно, проще всего минимизировать разности между фактическими значениями Yi, и предсказанными значениями . Однако, поскольку эти разности могут быть как положительными, так и отрицательными, следует минимизировать сумму их квадратов.

Поскольку = b0 + b1Xi, сумма квадратов принимает следующий вид:

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

До того, как Excel взял на себя всю рутинную работу, вычисления по методу наименьших квадратов были очень трудоемкими. Excel позволяет решать подобные задачи двумя способами. Во-первых, можно воспользоваться Пакетом анализа (строка Регрессия). Результаты представлены на рис. 4. Во-вторых, можно, выделив точки на графике (как на рис. 3б), кликнуть правой кнопкой мыши и выбрать Добавить линию тренда. Далее можно выбрать вид линии тренда (в нашем случае – Линейная), отформатировать линию, показать на графике уравнение и величину достоверности аппроксимации (R2) (рис. 5).

Рис. 4. Результаты решения задачи о зависимости между площадями и годовыми объемами продаж в магазинах сети Sunflower (получены с помощью Пакета анализа Excel)

Рис. 5. Диаграмма разброса и линия регрессии (тренда) в задаче о выборе магазина

Как следует из рис. 4 и 5, b0 = 0,9645, а b1 = 1,6699. Таким образом, уравнение линейной регрессии для этих данных имеет следующий вид: = 0,9645 + 1,6699Xi. Вычисленный наклон b1 = +1,6699. Это означает, что при возрастании переменной X на единицу среднее значение переменной Y возрастает на 1,6699 единиц. Иначе говоря, увеличение площади магазина на один квадратный фут приводит к увеличению годового объема продаж на 1,67 тыс. долл. Таким образом, наклон представляет собой долю годового объема продаж, зависящую от размера магазина. Вычисленный сдвиг b0 = +0,9645 (млн. долл.). Эта величина представляет собой среднее значение переменной Y при X = 0. Поскольку площадь магазина не может равняться нулю, сдвиг можно считать долей годового дохода, зависящей от других факторов. Следует отметить, однако, что сдвиг переменной Y выходит за пределы диапазона переменной X. Следовательно, к интерпретации параметра b0 необходимо относиться внимательно.

Пример 1. Один экономист решил предсказать изменение индекса 500 наиболее активно покупаемых акций на Нью-Йоркской фондовой бирже, публикуемого агентством Standard and Poor, на основе показателей экономики США за 50 лет. В результате он получил следующее уравнение линейной регрессии: Ŷi = –5,0 + 7Хi. Какой смысл имеют параметры сдвига b0 и наклона b1.

Решение. Сдвиг регрессии b0 равен –5,0. Это значит, что если рост экономики США равен нулю, индекс акций за год снизится на 5%. Наклон b1 равен 7. Следовательно, при увеличении темпов роста экономики на 1% индекс акций возрастает на 7%.

Пример 2. Вернемся к сценарию, изложенному в начале заметки. Применим модель линейной регрессии для прогноза объема годовых продаж во всех новых магазинах в зависимости от их размеров. Предположим, что площадь магазина равна 4000 квадратных футов. Какой среднегодовой объем продаж можно прогнозировать?

Решение. Подставим значение X = 4 (тыс. кв. футов) в уравнение линейной регрессии: = 0,9645 + 1,6699Xi = 0,9645 + 1,6699*4 = 7,644 млн. долл. Итак, прогнозируемый среднегодовой объем продаж в магазине, площадь которого равна 4000 кв. футов, составляет 7 644 000 долл.

Прогнозирование в регрессионном анализе: интерполяция и экстраполяция

Применяя регрессионную модель для прогнозирования, необходимо учитывать лишь допустимые значения независимой переменной. В этот диапазон входят все значения переменной X, начиная с минимальной и заканчивая максимальной. Таким образом, предсказывая значение переменной Y при конкретном значении переменной X, исследователь выполняет интерполяцию между значениями переменной X в диапазоне возможных значений. Однако экстраполяция значений за пределы этого интервала не всегда релевантна. Например, пытаясь предсказать среднегодовой объем продаж в магазине, зная его площадь (рис. 3а), мы можем вычислять значение переменной Y лишь для значений X от 1,1 до 5,8 тыс. кв. футов. Следовательно, прогнозировать среднегодовой объем продаж можно лишь для магазинов, площадь которых не выходит за пределы указанного диапазона. Любая попытка экстраполяции означает, что мы предполагаем, будто линейная регрессия сохраняет свой характер за пределами допустимого диапазона.

Оценки изменчивости

Вычисление сумм квадратов. Для того чтобы предсказать значение зависимой переменной по значениям независимой переменной в рамках избранной статистической модели, необходимо оценить изменчивость. Существует несколько способов оценки изменчивости. Первый способ использует общую сумму квадратов (total sum of squares — SST), позволяющую оценить колебания значений Yi вокруг среднего значения . В регрессионном анализе полная вариация, представляющая собой полную сумму квадратов, разделяется на объяснимую вариацию, или сумму квадратов регрессии (regression sum of squares — SSR), и необъяснимую вариацию, или сумму квадратов ошибок (error sum of squares — SSE). Объяснимая вариация характеризует взаимосвязь между переменными X и Y, а необъяснимая зависит от других факторов (рис. 6).

Рис. 6. Оценки изменчивости в модели регрессии

Сумма квадратов регрессии (SSR) представляет собой сумму квадратов разностей между Ŷi (предсказанным значением переменной Y) и (средним значением переменной Y). Сумма квадратов ошибок (SSE) является частью вариации переменной Y, которую невозможно описать с помощью регрессионной модели. Эта величина зависит от разностей между наблюдаемыми и предсказанными значениями.

Полная сумма квадратов (SST) равна сумме квадратов регрессии плюс сумма квадратов ошибок:

(3)   SST = SSR + SSE

Полная сумма квадратов (SST) равна сумме квадратов разностей между наблюдаемыми значениями переменной Y и ее средним значением:

Сумма квадратов регрессии (SSR) равна сумме квадратов разностей между предсказанными значениями переменной Y и ее средним значением:

Сумма квадратов ошибок (SSE) равна сумме квадратов разностей между наблюдаемыми и предсказанными значениями переменной Y:

Суммы квадратов, вычисленные с помощью программы Пакета анализа Excel при решении задачи о сети магазинов Sunflowers, представлены на рис. 4.

Полная сумма квадратов разностей равна SST = 116,9543. Эта величина состоит из суммы квадратов регрессии (SSR) равной 105,7476, и суммы квадратов ошибок (SSE), равной 11,2067.

Коэффициент смешанной корреляции. Величины SSR, SSE и SST не имеют очевидной интерпретации. Однако отношение суммы квадратов регрессии (SSR) к полной сумме квадратов (SST) представляет собой оценку полезности регрессионного уравнения. Это отношение называется коэффициентом смешанной корреляции r2:

Коэффициент смешанной корреляции оценивает долю вариации переменной Y, которая объясняется независимой переменной X в регрессионной модели. В задаче о сети магазинов Sunflowers SSR = 105,7476 и SST = 116,9543. Следовательно, r2 = 105,7476 / 116,9543 = 0,904. Таким образом, 90,4% вариации годового объема продаж объясняется изменчивостью площади магазинов, измеренной в квадратных футах. Данная величина r2 свидетельствует о сильной положительной линейной взаимосвязи между двумя переменными, поскольку применение регрессионной модели снижает изменчивость прогнозируемых годовых объемов продаж на 90,4%. Только 9,6% изменчивости годовых объемов продаж в выборке магазинов объясняются другими факторами, не учтенными в регрессионной модели.

Коэффициент смешанной корреляции в задаче о сети магазинов Sunflowers представлен в таблице Регрессионная статистика на рис. 4.

Среднеквадратичная ошибка оценки. Хотя метод наименьших квадратов позволяет вычислить линию, минимизирующую отклонение от наблюдаемых значений, наличие суммы квадратов ошибок (SSE) свидетельствует о том, что линейная регрессия не дает абсолютной точности прогноза, если, конечно, точки наблюдения не лежат на регрессионной прямой. Однако ожидать этого так же неестественно, как предполагать, что все выборочные значения точно равны их среднему арифметическому. Следовательно, необходима статистика, которая позволила бы оценить отклонение предсказанных значений переменной Y от ее реальных значений, аналогично тому, как стандартное отклонение, введенное ранее, позволяет оценить колебание данных вокруг их средней величины. Стандартное отклонение наблюдаемых значений переменной Y от ее регрессионной прямой называется среднеквадратичной ошибкой оценки. Отклонение реальных данных от регрессионной прямой в задаче о сети магазинов Sunflowers показано на рис. 5.

Среднеквадратичная ошибка оценки

где Yi — фактическое значение переменной Y при заданном значении Xi, Ŷi — предсказанное значение переменной Y при заданном значении Xi, SSE — сумма квадратов ошибок.

Поскольку SSE = 11,2067, по формуле (8) получаем:

Таким образом, среднеквадратичная ошибка оценки равна 0,9664 млн. долл. (т.е. 966 400 долл.). Этот параметр также рассчитывается Пакетом анализа (см. рис. 4). Среднеквадратичная ошибка оценки характеризует отклонение реальных данных от линии регрессии. Она измеряется в тех же единицах, что и переменная Y. По смыслу среднеквадратичная ошибка очень похожа на стандартное отклонение. В то время как стандартное отклонение характеризует разброс данных вокруг их среднего значения, среднеквадратичная ошибка позволяет оценить колебание точек наблюдения вокруг регрессионной прямой. Cреднеквадратичная ошибка оценки позволяет обнаружить статистически значимую зависимость, существующую между двумя переменными, и предсказать значения переменной Y.

Предположения

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

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

Первое предположение, о нормальном распределении ошибок, требует, чтобы при каждом значении переменной X ошибки линейной регрессии имели нормальное распределение (рис. 7). Как и t— и F-критерий дисперсионного анализа, регрессионный анализ довольно устойчив к нарушениям этого условия. Если распределение ошибок относительно линии регрессии при каждом значении X не слишком сильно отличается от нормального, выводы относительно линии регрессии и коэффициентов регрессии изменяются незначительно.

Рис. 7. Предположение о нормальном распределении ошибок

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

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

Анализ остатков

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

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

(9) ei = Yi – Ŷi

Для оценки пригодности эмпирической модели регрессии остатки откладываются по вертикальной оси, а значения Xi — по горизонтальной. Если эмпирическая модель пригодна, график не должен иметь ярко выраженной закономерности. Если же модель регрессии не пригодна, на рисунке проявится зависимость между значениями Xi и остатками еi.

Рассмотрим примеры (рис. 8). Панель А иллюстрирует возрастание переменной Y при увеличении переменной X. Однако зависимость между этими переменными носит нелинейный характер, поскольку скорость возрастания переменной Y падает при увеличении переменной X. Таким образом, для аппроксимации зависимости между этими переменными лучше подойдет квадратичная модель. Особенно ярко квадратичная зависимость между величинами Xi и ei проявляется на панели Б. Графическое изображение остатков позволяет отфильтровать или удалить линейную зависимость между переменными X и Y и выявить недостаточную точность модели простой линейной регрессии. Таким образом, в данной ситуации вместо простой линейной модели должна применяться квадратичная модель, обладающая более высокой точностью.

Рис. 8. Исследование эмпирической модели простой линейной регрессии

Вернемся к задаче о сети магазинов Sunflowers и посмотрим, хорошо ли подходит простая линейная регрессия для ее решения. Соответствующие данные и расчеты приведены на рис. 9а (формулы можно посмотреть в Excel-файле). Построим диаграмму разброса, откладывая по вертикальной оси остатки ei, а по горизонтальной — независимую переменную Xi (рис. 9б). Несмотря на большой разброс остатков, между ei и Хi нет ярко выраженной зависимости. Остатки одинаково часто принимают как положительные, так и отрицательные значения. Это позволяет сделать вывод, что модель линейной регрессии пригодна для решения задачи о сети магазинов Sunflowers.

Рис. 9. Остатки ei, вычисленные при решении задачи о сети магазинов Sunflowers

Значения остатков (таблица на рис. 9а) и график остатков (аналог рис. 9б) можно получить непосредственно в процедуре Регрессия Пакета анализа. Просто поставьте соответствующие галки (рис. 10).

Рис. 10. Остатки ei и график остатков полученные с помощью Пакета анализа

Проверка условий. График остатков позволяет оценить вариации ошибок. На рис. 10 нет особых различий между ошибками, соответствующими разным значениям Xi. Следовательно, вариации ошибок при разных значениях Хi приблизительно одинаковы. Рассмотрим гипотетическую ситуацию, в которой это условие не выполняется (рис. 11). На этом рисунке изображен эффект веера: при возрастании значений Хi ошибки увеличиваются. Таким образом, изменчивость значений Yi при разных значениях Хi является непостоянной.

Рис. 11. Пример нарушения условия независимости вариаций ошибок от Xi

Нормальность. Чтобы проверить предположение о нормальном распределении ошибок, построим график нормального распределения на основе точечного графика, на вертикальной оси которого отложены значения остатков, а на горизонтальной оси — соответствующие квантили стандартизованного нормального распределения (подробнее см. Проверка гипотезы о нормальном распределении). Для построения такого графика значения остатков должны быть упорядочены по возрастанию (рис. 12). График нормального распределения может быть построен одним кликом с помощью Пакета анализа Excel – просто поставьте соответствующую галочку в окне Регрессия (см. рис. 10, самый низ окна Регрессия – опция График нормальной вероятности).

Рис. 12. График нормального распределения для остатков

Без визуализации данных (с помощью гистограммы, диаграммы «ствол и листья», блочной диаграммы или графика как на рис. 12) проверить предположение о нормальном распределении ошибок очень трудно. Данные, изображенные на рис. 12, не слишком сильно отличаются от нормального распределения. Устойчивость регрессионного анализа и небольшой объем выборки позволяют утверждать, что условие о нормальном распределении ошибок нарушается незначительно.

Независимость. Предположение о независимости ошибок также проверяется с помощью графика остатков. Данные, собранные на протяжении некоторого периода времени, иногда демонстрируют эффект автокорреляции между последовательными наблюдениями. В таких ситуациях остатки зависят от значений предыдущих остатков. Подобная связь между остатками нарушает предположение о независимости ошибок. Эффект автокорреляции хорошо выявляется на графике. Кроме того, его можно измерить с помощью процедуры Дурбина-Уотсона (см. ниже). Если данные о размерах магазинов и объемах продаж собирались в течение одного и того же периода времени, гипотезу об их независимости проверять не имеет смысла.

Измерение автокорреляции: статистика Дурбина–Уотсона

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

Распознавание автокорреляции с помощью графика остатков. Для выявления автокорреляции необходимо упорядочить остатки по времени и построить их график. Если данные обладают положительной автокорреляцией, на графике возникнут кластеры остатков, имеющие одинаковый знак. В случае отрицательной автокорреляции остатки будут скачкообразно принимать то положительные, то отрицательные значения. Этот вид автокорреляции очень редко встречается в регрессионном анализе, поэтому мы рассмотрим лишь положительную автокорреляцию. Проиллюстрируем ее следующим примером. Предположим, что менеджер магазина, доставляющего товары на дом, пытается предсказать объем продаж по количеству клиентов, совершивших покупки в течение 15 недель (рис. 13).

Рис. 13. Количество клиентов и объемы продаж за 15 недель

Поскольку данные собирались на протяжении 15 последовательных недель в одном и том же магазине, необходимо определить, наблюдается ли эффект автокорреляции. Построим регрессию с использованием Пакета анализа; включим вывод Остатков, но не будем включать График остатков (рис. 14).

Рис. 14. Параметры линейной регрессии, полученные с использованием Пакета анализа

Анализ рис. 14 показывает, что r2 = 0,657. Это значит, что 65,7% вариации объемов продаж объясняется изменчивостью количества клиентов. Кроме того, сдвиг b0 переменной Y равен –16,032, а наклон b1 = 0,0308. Однако, прежде чем применять эту модель, необходимо выполнить анализ остатков. Поскольку данные собирались на протяжении 15 последовательных недель, их следует отобразить на графике в том же порядке (рис. 15).

Рис. 15. Зависимость остатков от времени

Анализ рис. 15 показывает, что остатки циклически колеблются вверх и вниз. Эта цикличность является явным признаком автокорреляции. Следовательно, гипотезу о независимости остатков следует отклонить.

Статистика Дурбина-Уотсона. Автокорреляцию можно выявить и измерить с помощью статистики Дурбина-Уотсона. Эта статистика оценивает корреляцию между соседними остатками:

где еi — остаток, соответствующий i-му периоду времени.

Чтобы лучше понять статистику Дурбина-Уотсона, рассмотрим ее составные части. Числитель представляет собой сумму квадратов разностей между соседними остатками, начиная со второго и заканчивая n-м наблюдением. Знаменатель является суммой квадратов остатков. Вот, что по этому поводу написано в Википедии:

где ρ1 – коэффициент автокорреляции; если ρ1 = 0 (нет автокорреляции), D ≈ 2; если ρ1 ≈ 1 (положительная автокорреляции), D ≈ 0; если ρ1 = -1 (отрицательная автокорреляции), D ≈ 4.

На практике применение критерия Дурбина-Уотсона основано на сравнении величины D с критическими теоретическими значениями dL и dU для заданного числа наблюдений n, числа независимых переменных модели k (для простой линейной регрессии k = 1) и уровня значимости α. Если D < dL, гипотеза о независимости случайных отклонений отвергается (следовательно, присутствует положительная автокорреляция); если D > dU, гипотеза не отвергается (то есть автокорреляция отсутствует); если dL < D < dU, нет достаточных оснований для принятия решения. Когда расчётное значение D превышает 2, то с dL и dU сравнивается не сам коэффициент D, а выражение (4 – D).

Для вычисления статистики Дурбина-Уотсона в Excel обратимся к нижней таблице на рис. 14 Вывод остатка. Числитель в выражении (10) вычисляется с помощью функции =СУММКВРАЗН(массив1;массив2), а знаменатель =СУММКВ(массив) (рис. 16).

Рис. 16. Формулы расчета статистики Дурбина-Уотсона

В нашем примере D = 0,883. Основной вопрос заключается в следующем — какое значение статистики Дурбина-Уотсона следует считать достаточно малым, чтобы сделать вывод о существовании положительной автокорреляции? Необходимо соотнести значение D с критическими значениями (dL и dU), зависящими от числа наблюдений n и уровня значимости α (рис. 17).

Рис. 17. Критические значения статистики Дурбина-Уотсона (фрагмент таблицы)

Таким образом, в задаче об объеме продаж в магазине, доставляющем товары на дом, существуют одна независимая переменная (k = 1), 15 наблюдений (n = 15) и уровень значимости α = 0,05. Следовательно, dL= 1,08 и dU = 1,36. Поскольку D = 0,883 < dL= 1,08, между остатками существует положительная автокорреляция, метод наименьших квадратов применять нельзя.

Проверка гипотез о наклоне и коэффициенте корреляции

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

Применение t-критерия для наклона. Проверяя, равен ли наклон генеральной совокупности β1 нулю, можно определить, существует ли статистически значимая зависимость между переменными X и Y. Если эта гипотеза отклоняется, можно утверждать, что между переменными X и Y существует линейная зависимость. Нулевая и альтернативная гипотезы формулируются следующим образом: Н0: β1 = 0 (нет линейной зависимости), Н1: β1 ≠ 0 (есть линейная зависимость). По определению t-статистика равна разности между выборочным наклоном и гипотетическим значением наклона генеральной совокупности, деленной на среднеквадратичную ошибку оценки наклона:

(11) t = (b1β1) / Sb1

где b1 – наклон прямой регрессии по выборочным данным, β1 – гипотетический наклон прямой генеральной совокупности, , а тестовая статистика t имеет t-распределение с n – 2 степенями свободы.

Проверим, существует ли статистически значимая зависимость между размером магазина и годовым объемом продаж при α = 0,05. t-критерий выводится наряду с другими параметрами при использовании Пакета анализа (опция Регрессия). Полностью результаты работы Пакета анализа приведены на рис. 4, фрагмент, относящийся к t-статистике – на рис. 18.

Рис. 18. Результаты применения t-критерия, полученные с помощью Пакета анализа Excel

Поскольку число магазинов n = 14 (см. рис.3), критическое значение t-статистики при уровне значимости α = 0,05 можно найти по формуле: tL =СТЬЮДЕНТ.ОБР(0,025;12) = –2,1788, где 0,025 – половина уровня значимости, а 12 = n – 2; tU =СТЬЮДЕНТ.ОБР(0,975;12) = +2,1788.

Поскольку t-статистика = 10,64 > tU = 2,1788 (рис. 19), нулевая гипотеза Н0 отклоняется. С другой стороны, р-значение для Х = 10,6411, вычисляемое по формуле =1-СТЬЮДЕНТ.РАСП(D3;12;ИСТИНА), приближенно равно нулю, поэтому гипотеза Н0 снова отклоняется. Тот факт, что р-значение почти равно нулю, означает, что если бы между размерами магазинов и годовым объемом продаж не существовало реальной линейной зависимости, обнаружить ее с помощью линейной регрессии было бы практически невозможно. Следовательно, между средним годовым объемом продаж в магазинах и их размером существует статистически значимая линейная зависимость.

Рис. 19. Проверка гипотезы о наклоне генеральной совокупности при уровне значимости, равном 0,05, и 12 степенях свободы

Применение F-критерия для наклона. Альтернативным подходом к проверке гипотез о наклоне простой линейной регрессии является использование F-критерия. Напомним, что F-критерий применяется для проверки отношения между двумя дисперсиями (подробнее см. Однофакторный дисперсионный анализ). При проверке гипотезы о наклоне мерой случайных ошибок является дисперсия ошибки (сумма квадратов ошибок, деленная на количество степеней свободы), поэтому F-критерий использует отношение дисперсии, объясняемой регрессией (т.е. величины SSR, деленной на количество независимых переменных k), к дисперсии ошибок (MSE = SYX2).

По определению F-статистика равна среднему квадрату отклонений, обусловленных регрессией (MSR), деленному на дисперсию ошибки (MSE): F = MSR/MSE, где MSR = SSR / k, MSE = SSE/(n– k – 1), k – количество независимых переменных в регрессионной модели. Тестовая статистика F имеет F-распределение с k и n – k – 1 степенями свободы.

При заданном уровне значимости α решающее правило формулируется так: если F > FU, нулевая гипотеза отклоняется; в противном случае она не отклоняется. Результаты, оформленные в виде сводной таблицы дисперсионного анализа, приведены на рис. 20.

Рис. 20. Таблица дисперсионного анализа для проверки гипотезы о статистической значимости коэффициента регрессии

Аналогично t-критерию F-критерий выводится в таблицу при использовании Пакета анализа (опция Регрессия). Полностью результаты работы Пакета анализа приведены на рис. 4, фрагмент, относящийся к F-статистике – на рис. 21.

Рис. 21. Результаты применения F-критерия, полученные с помощью Пакета анализа Excel

F-статистика равна 113,23, а р-значение близко к нулю (ячейка Значимость F). Если уровень значимости α равен 0,05, определить критическое значение F-распределения с одной и 12 степенями свободы можно по формуле FU =F.ОБР(1-0,05;1;12) = 4,7472 (рис. 22). Поскольку F = 113,23 > FU = 4,7472, причем р-значение близко к 0 < 0,05, нулевая гипотеза Н0 отклоняется, т.е. размер магазина тесно связан с его годовым объемом продаж.

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

Доверительный интервал, содержащий наклон β1. Для проверки гипотезы о существовании линейной зависимости между переменными можно построить доверительный интервал, содержащий наклон β1 и убедиться, что гипотетическое значение β1 = 0 принадлежит этому интервалу. Центром доверительного интервала, содержащего наклон β1, является выборочный наклон b1, а его границами — величины b1 ± tn–2Sb1

Как показано на рис. 18, b1 = +1,670, n = 14, Sb1 = 0,157. t12 =СТЬЮДЕНТ.ОБР(0,975;12) = 2,1788. Следовательно, b1 ± tn–2Sb1 = +1,670 ± 2,1788 * 0,157 = +1,670 ± 0,342, или + 1,328 ≤ β1 ≤ +2,012. Таким образом, наклон генеральной совокупности с вероятностью 0,95 лежит в интервале от +1,328 до +2,012 (т.е. от 1 328 000 до 2 012 000 долл.). Поскольку эти величины больше нуля, между годовым объемом продаж и площадью магазина существует статистически значимая линейная зависимость. Если бы доверительный интервал содержал нуль, между переменными не было бы зависимости. Кроме того, доверительный интервал означает, что каждое увеличение площади магазина на 1 000 кв. футов приводит к увеличению среднего объема продаж на величину от 1 328 000 до 2 012 000 долларов.

Использование t-критерия для коэффициента корреляции. Ранее был введен коэффициент корреляции r, представляющий собой меру зависимости между двумя числовыми переменными. С его помощью можно установить, существует ли между двумя переменными статистически значимая связь. Обозначим коэффициент корреляции между генеральными совокупностями обеих переменных символом ρ. Нулевая и альтернативная гипотезы формулируются следующим образом: Н0: ρ = 0 (нет корреляции), Н1: ρ ≠ 0 (есть корреляция). Проверка существования корреляции:

где r = +, если b1 > 0, r = –, если b1 < 0. Тестовая статистика t имеет t-распределение с n – 2 степенями свободы.

В задаче о сети магазинов Sunflowers r2 = 0,904, а b1— +1,670 (см. рис. 4). Поскольку b1 > 0, коэффициент корреляции между объемом годовых продаж и размером магазина равен r = +√0,904 = +0,951. Проверим нулевую гипотезу, утверждающую, что между этими переменными нет корреляции, используя t-статистику:

При уровне значимости α = 0,05 нулевую гипотезу следует отклонить, поскольку t = 10,64 > 2,1788. Таким образом, можно утверждать, что между объемом годовых продаж и размером магазина существует статистически значимая связь.

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

Оценка математического ожидания и предсказание индивидуальных значений

В этом разделе рассматриваются методы оценки математического ожидания отклика Y и предсказания индивидуальных значений Y при заданных значениях переменной X.

Построение доверительного интервала. В примере 2 (см. выше раздел Метод наименьших квадратов) регрессионное уравнение позволило предсказать значение переменной Y при заданном значении переменной X. В задаче о выборе места для торговой точки средний годовой объем продаж в магазине площадью 4000 кв. футов был равен 7,644 млн. долл. Однако эта оценка математического ожидания генеральной совокупности является точечной. Ранее для оценки математического ожидания генеральной совокупности была предложена концепция доверительного интервала. Аналогично можно ввести понятие доверительного интервала для математического ожидания отклика при заданном значении переменной X:

где , = b0 + b1Xi – предсказанное значение переменное Y при X = Xi, SYX – среднеквадратичная ошибка, n – объем выборки, Xi — заданное значение переменной X, µY|X=Xi – математическое ожидание переменной Y при Х = Хi, SSX =

Анализ формулы (13) показывает, что ширина доверительного интервала зависит от нескольких факторов. При заданном уровне значимости возрастание амплитуды колебаний вокруг линии регрессии, измеренное с помощью среднеквадратичной ошибки, приводит к увеличению ширины интервала. С другой стороны, как и следовало ожидать, увеличение объема выборки сопровождается сужением интервала. Кроме того, ширина интервала изменяется в зависимости от значений Xi. Если значение переменной Y предсказывается для величин X, близких к среднему значению , доверительный интервал оказывается уже, чем при прогнозировании отклика для значений, далеких от среднего.

Допустим, что, выбирая место для магазина, мы хотим построить 95%-ный доверительный интервал для среднего годового объема продаж во всех магазинах, площадь которых равна 4000 кв. футов:

Следовательно, средний годовой объем продаж во всех магазинах, площадь которых равна 4 000 кв. футов, с 95% -ной вероятностью лежит в интервале от 6,971 до 8,317 млн. долл.

Вычисление доверительного интервала для предсказанного значения. Кроме доверительного интервала для математического ожидания отклика при заданном значении переменной X, часто необходимо знать доверительный интервал для предсказанного значения. Несмотря на то что формула для вычисления такого доверительного интервала очень похожа на формулу (13), этот интервал содержит предсказанное значение, а не оценку параметра. Интервал для предсказанного отклика YX=Xi при конкретном значении переменной Xi определяется по формуле:

Предположим, что, выбирая место для торговой точки, мы хотим построить 95%-ный доверительный интервал для предсказанного годового объема продаж в магазине, площадь которого равна 4000 кв. футов:

Следовательно, предсказанный годовой объем продаж в магазине, площадь которого равна 4000 кв. футов, с 95%-ной вероятностью лежит в интервале от 5,433 до 9,854 млн. долл. Как видим, доверительный интервал для предсказанного значения отклика намного шире, чем доверительный интервал для его математического ожидания. Это объясняется тем, что изменчивость при прогнозировании индивидуальных значений намного больше, чем при оценке математического ожидания.

Подводные камни и этические проблемы, связанные с применением регрессии

Трудности, связанные с регрессионным анализом:

  • Игнорирование условий применимости метода наименьших квадратов.
  • Ошибочная оценка условий применимости метода наименьших квадратов.
  • Неправильный выбор альтернативных методов при нарушении условий применимости метода наименьших квадратов.
  • Применение регрессионного анализа без глубоких знаний о предмете исследования.
  • Экстраполяция регрессии за пределы диапазона изменения объясняющей переменной.
  • Путаница между статистической и причинно-следственной зависимостями.

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

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

Рис. 23. Четыре набора искусственных данных

Рис. 24. Регрессионный анализ четырех искусственных наборов данных; выполнен с помощью Пакета анализа (кликните на рисунке, чтобы увеличить изображение)

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

Рис. 25. Диаграммы разброса для четырех наборов данных

Диаграммы разброса и графики остатков свидетельствуют о том, что эти данные отличаются друг от друга. Единственный набор, распределенный вдоль прямой линии, — набор А. График остатков, вычисленных по набору А, не имеет никакой закономерности. Этого нельзя сказать о наборах Б, В и Г. График разброса, построенный по набору Б, демонстрирует ярко выраженную квадратичную модель. Этот вывод подтверждается графиком остатков, имеющим параболическую форму. Диаграмма разброса и график остатков показывают, что набор данных В содержит выброс. В этой ситуации необходимо исключить выброс из набора данных и повторить анализ. Метод, позволяющий обнаруживать и исключать выбросы из наблюдений, называется анализом влияния. После исключения выброса результат повторной оценки модели может оказаться совершенно иным. Диаграмма разброса, построенная по данным из набора Г, иллюстрирует необычную ситуацию, в которой эмпирическая модель значительно зависит от отдельного отклика (Х8 = 19, Y8 = 12,5). Такие регрессионные модели необходимо вычислять особенно тщательно. Итак, графики разброса и остатков являются крайне необходимым инструментом регрессионного анализа и должны быть его неотъемлемой частью. Без них регрессионный анализ не заслуживает доверия.

Рис. 26. Графики остатков для четырех наборов данных

Как избежать подводных камней при регрессионном анализе:

  • Анализ возможной взаимосвязи между переменными X и Y всегда начинайте с построения диаграммы разброса.
  • Прежде чем интерпретировать результаты регрессионного анализа, проверяйте условия его применимости.
  • Постройте график зависимости остатков от независимой переменной. Это позволит определить, насколько эмпирическая модель соответствует результатам наблюдения, и обнаружить нарушение постоянства дисперсии.
  • Для проверки предположения о нормальном распределении ошибок используйте гистограммы, диаграммы «ствол и листья», блочные диаграммы и графики нормального распределения.
  • Если условия применимости метода наименьших квадратов не выполняются, используйте альтернативные методы (например, модели квадратичной или множественной регрессии).
  • Если условия применимости метода наименьших квадратов выполняются, необходимо проверить гипотезу о статистической значимости коэффициентов регрессии и построить доверительные интервалы, содержащие математическое ожидание и предсказанное значение отклика.
  • Избегайте предсказывать значения зависимой переменной за пределами диапазона изменения независимой переменной.
  • Имейте в виду, что статистические зависимости не всегда являются причинно-следственными. Помните, что корреляция между переменными не означает наличия причинно-следственной зависимости между ними.

Резюме. Как показано на структурной схеме (рис. 27), в заметке описаны модель простой линейной регрессии, условия ее применимости и способы проверки этих условий. Рассмотрен t-критерий для проверки статистической значимости наклона регрессии. Для предсказания значений зависимой переменной использована регрессионная модель. Рассмотрен пример, связанный с выбором места для торговой точки, в котором исследуется зависимость годового объема продаж от площади магазина. Полученная информация позволяет точнее выбрать место для магазина и предсказать его годовой объем продаж. В следующих заметках будет продолжено обсуждение регрессионного анализа, а также рассмотрены модели множественной регрессии.

Рис. 27. Структурная схема заметки

Предыдущая заметка Критерий согласия «хи-квадрат»

Следующая заметка Введение в множественную регрессию

К оглавлению Статистика для менеджеров с использованием Microsoft Excel


[1] Используются материалы книги Левин и др. Статистика для менеджеров. – М.: Вильямс, 2004. – с. 792–872

[2] Если зависимая переменная является категорийной, необходимо применять логистическую регрессию.

Регрессия как задача машинного обучения

38 мин на чтение

(55.116 символов)

Постановка задачи регрессии

Задача регрессии
Источник: Analytics Vidhya.

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

Задача регрессии относится к категории задач обучения с учителем. Это значит, что набор данных, который используется для обучения, должен иметь определенную структуру. Обычно, наборы данных для машинного обучения представляют собой таблицу, в которой по строкам перечислены разные объекты наблюдений или измерений. В столбцах — различные характеристики, или атрибуты, объектов. А на пересечении строк и столбцов — значение данной характеристики у данного объекта. Обычно один атрибут (или переменная) имеет особый характер — именно ее значение мы и хотим научиться предсказывать с помощью модели машинного обучения. Эта характеристика объекта называется целевая переменная. И если эта целевая переменная выражена числом (а точнее, некоторой непрерывной величиной) — то мы говорим о задаче регрессии.

Задачи регрессии на практике встречаются довольно часто. Например, предсказание цены объекта недвижимости — классическая регрессионная задача. В таких проблемах атрибутами выступают разные характеристики квартир или домов — площадь, этажность, расположение, расстояние до центра города, количество комнат, год постройки. В разных наборах данных собрана разная информация И, соответственно, модели тоже должны быть разные. Другой пример — предсказание цены акций или других финансовых активов. Или предсказание температуры завтрашним днем.

Во всех таких задачах нам нужно иметь данные, которые позволят осуществить такое предсказание. Да, “предсказание” — это условный термин, не всегда мы говорим о будущих событиях. Регрессионные модели используют информацию об объектах в обучающем наборе данных, чтобы сделать вывод о возможном значении целевой переменной. И для этого нужно, чтобы ее значение имело какую-то зависимость от имеющихся у нас атрибутов. Если построить модель предсказания цены акции, но на вход подать информацию о футбольных матчах — ничего не получится. Мы предполагаем, что в наборе данных собраны именно те атрибуты объектов, которые имеют влияние на на значение целевой переменной. И чем больше это предположение выполняется, тем точнее будет потенциально наша модель.

Немного поговорим о терминах. Набор данных который мы используем для обучения модели называют датасетом (dataset) или обучающей выборкой (training set). Объекты, которые описываются в датасете еще называют точками данных (data points). Целевую переменную еще называют на статистический манер зависимой переменной (dependent variable) или результативной, выходной (output), а остальные атрибуты — независимыми переменными (dependent variables), или признаками (features), или факторами, или входными переменными (input). Значения одного конкретного атрибута для всех объектов обучающей выборки часто представляют как вектор этого признака (feature vector). А всю таблицу всех атрибутов называют матрицей атрибутов (feature matrix). Соответственно, еще есть вектор целевой переменной, он не входит в матрицу атрибутов.

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

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

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

Выводы:

  1. Регрессия — это задача машинного обучения с учителем, которая заключается в предсказании некоторой непрерывной величины.
  2. Для использования регрессионных моделей нужно, чтобы в датасете были характеристики объектов и “правильные” значения целевой переменной.
  3. Примеры регрессионных задач — предсказание цены акции, оценка цены объекта недвижимости.
  4. Задача регрессии основывается на предположении, что значение целевой переменной зависит от значения признаков.
  5. Регрессионная модель принимает набор значений и выдает предсказание значения целевой переменной.
  6. В качестве регрессионных моделей часто берут аналитические функции, например, линейную.

Линейная регрессия с одной переменной

Функция гипотезы

Модель регрессии

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

Парная, или одномерная (univariate) регрессия используется, когда вы хотите предсказать одно выходное значение (чаще всего обозначаемое $y$), зависящее от одного входного значения (обычно обозначается $x$). Сама функция называется функцией гипотезы или моделью. В качестве функции гипотезы для парной регрессии можно выбрать любую функцию, но мы пока потренируемся с самой простой функцией одной переменной — линейной функцией. Тогда нашу модель можно назвать парной линейной регрессией.

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

[hat{y} = h_b (x) = b_0 + b_1 x]

Обратите внимание, что это похоже на уравнение прямой. Эта модель соответствует множеству всех возможных прямых на плоскости. Когда мы конкретизируем модель значениями параметров (в данном случае — $b_0$ и $b_1$), мы получаем конкретную прямую. И наша задача состоит в том, чтобы выбрать такую прямую, которая бы лучше всего “легла” в точки из нашей обучающей выборки.

В данном случае, мы пытаемся подобрать функцию h(x) таким образом, чтобы отобразить данные нам значения x в данные значения y.

Допустим, мы имеем следующий обучающий набор данных:

входная переменная x выходная переменная y
0 4
1 7
2 7
3 8

Мы можем составить случайную гипотезу с параметрами $ b_0 = 2, b_1 = 2 $. Тогда для входного значения $ x=1 $ модель выдаст предсказание, что $ y=4 $, что на 3 меньше данного. Значение $y$б которое посчитала модель будем называть теоретическим или предсказанным (predicted), а значение, которое дано в наборе данных — эмпирическим или истинным (true). Задача регрессии состоит в нахождении таких параметров функции гипотезы, чтобы она отображала входные значения в выходные как можно более точно, или, другими словами, описывала линию, наиболее точно ложащуюся в данные точки на плоскости $(x, y)$.

Выводы:

  1. Модель машинного обучения — это параметрическая функция.
  2. Задача обучения состоит в том, чтобы подобрать параметры модели таким образом, чтобы она лучше всего описывала обучающие данные.
  3. Парная линейная регрессия работает, если есть всего одна входящая переменная.
  4. Парная линейная регрессия — одна из самых простых моделей машинного обучения.
  5. Парная линейная регрессия соответствует множеству всех прямых на плоскости. Из них мы выбираем одну, наиболее подходящую.

Функция ошибки

Как мы уже говорили, разные значения параметров дают разные модели. Для того, чтобы подобрать наилучшую модель, нам нужно средство измерения “точности” модели, некоторая функция, которая показывает, насколько модель хорошо или плохо соответствует имеющимся данным.

Разные модели

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

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

Отклонения значений

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

[J(b_0, b_1)
= frac{1}{2m} sum_{i=1}^{m} (hat{y_i} — y_i)^2
= frac{1}{2m} sum_{i=1}^{m} (h_b(x_i) — y_i)^2]

Эту функцию называют «функцией квадрата ошибки» или «среднеквадратичной ошибкой» (mean squared error, MSE). Среднее значение уменьшено вдвое для удобства вычисления градиентного спуска, так как производная квадратичной функции будет отменять множитель 1/2. Вообще, функцию ошибки можно свободно домножить или разделить на любое число (положительное), ведь нам не важна конкретная величина этой функции. Нам важно, что какие-то модели (то есть наборы значений параметров модели) имеют низкую ошибку, они нам подходят больше, а какие-то — высокую ошибку, они подходят нам меньше.

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

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

Давайте проследим формирование функции ошибки на еще более простом примере. Возьмем упрощенную форму линейной модели — прямую пропорциональность. Она выражается формулой:

[hat{y} = h_b (x) = b_1 x]

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

Функция ошибки одной переменной

При значении $b_1 = -1$ линия существенно отклоняется от точек. Отметим уровень ошибки (примерно 10) на правом графике.

Функция ошибки одной переменной

Если взять значение $b_1 = 0$ линия гораздо ближе к точкам, но ошибка все еще есть. Отметим новое значение на правом графике в точке 0.

Функция ошибки одной переменной

При значении $b_1 = 1$ график точно ложится в точки, таким образом ошибка становится равной нулю. Отмечаем ее так же.

Функция ошибки одной переменной

При дальнейшем увеличении $b_1$ линия становится выше точек. Но функция ошибки все равно будет положительной. Теперь она опять станет расти.

Функция ошибки одной переменной

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

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

Функция ошибки одной переменной

Функция ошибки одной переменной

Функция ошибки одной переменной

Функция ошибки одной переменной

Функция ошибки одной переменной

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

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

Среднеквадратическая ошибка

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

Если мы попытаемся представить это наглядно, наш набор данных обучения будет разбросан по плоскости x-y. Мы пытаемся подобрать прямую линию, которая проходит через этот разбросанный набор данных. Наша цель — получить наилучшую возможную линию. Лучшая линия будет такой, чтобы средние квадраты вертикальных расстояний точек от линии были наименьшими. В лучшем случае линия должна проходить через все точки нашего набора данных обучения. В таком случае значение J будет равно 0.

Ошибка

Ошибка

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

Выводы:

  1. Функция ошибки нужна для того, чтобы отличать хорошие модели от плохих.
  2. Функция ошибки показывает численно, насколько модель хорошо описывает данные.
  3. Аргументами функции ошибки являются параметры модели, ошибка зависит от них.
  4. Само значение функции ошибки не несет никакого смысла, оно используется только в сравнении.
  5. Цель алгоритма машинного обучения — минимизировать функцию ошибки, то есть найти такой набор параметров модели, при которых ошибка минимальна.
  6. Чаще всего используется так называемая L2-ошибка — средний квадрат отклонений теоретических значений от эмпирических (метрика MSE).

Метод градиентного спуска

Таким образом, у нас есть функция гипотезы, и способ оценить, насколько хорошо конкретная гипотеза вписывается в данные. Теперь нам нужно подобрать параметры функции гипотезы. Вот где приходит на помощь метод градиентного спуска.

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

Давайте найдем производную среднеквадратической функции ошибки:

[J(b_0, b_1) = frac{1}{2m} sum_{i=1}^{m} (h_b(x_i) — y_i)^2]

[J(b_0, b_1) = frac{1}{2m} sum_{i=1}^{m} (h_b(x_i) — y_i)^2]

[frac{partial}{partial b_i} J =
frac{1}{m} sum_{i=1}^{m} (h_b(x_i) — y^{(i)}) cdot frac{partial}{partial b_i} h_b(x_i)]

[J(b_0, b_1) = frac{1}{2m} sum_{i=1}^{m} (b_0 + b_1 x_i — y_i)^2]

[frac{partial J}{partial b_0} =
frac{1}{m} sum (b_0 + b_1 x_i — y_i) =
frac{1}{m} sum (h_b(x_i) — y_i)]

[frac{partial J}{partial b_1} =
frac{1}{m} sum (b_0 + b_1 x_i — y_i) cdot x_i =
frac{1}{m} sum (h_b(x_i) — y_i) cdot x_i]

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

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

Градиентный спуск

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

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

Алгоритм градиентного спуска:

повторяйте до сходимости:

[b_j := b_j — alpha frac{partial}{partial b_j} J(b_0, b_1)]

где j=0,1 — представляет собой индекс номера признака.

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

Алгоритм градиентного спуска имеет один параметр — скорость обучения. Он влияет на то, как быстро мы будем приближаться к оптимуму. Кажется, что чем быстрее, тем лучше, но оказывается, что если значение данного параметра слишком велико, то мы буем постоянно промахиваться и алгоритм будет расходиться.

Градиентный спуск

Алгоритм градиентного спуска для парной линейной регрессии:

повторяйте до сходимости:

[b_0 := b_0 — alpha frac{1}{m} sum_{i=1}^{m} (h_b(x^{(i)} )- y^{(i)})]

[b_1 := b_1 — alpha frac{1}{m} sum_{i=1}^{m} (h_b(x^{(i)}) — y^{(i)}) cdot x^{(i)}]

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

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

Спуск

Другой спуск

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

Выводы:

  1. Метод градиентного спуска нужен, чтобы найти минимум функции, если мы не можем ее вычислить аналитически.
  2. Это численный итеративный алгоритм локальной оптимизации.
  3. Для запуска градиентного спуска нужно знать частную производную функции ошибки.
  4. Для начала мы берем произвольные значения параметров, затем обновляем их по данной формуле.
  5. Доказано, что этот метод сходится к локальному минимуму.
  6. Если функция ошибки достаточно сложная, то разные начальные точки дадут разный результат.
  7. Метод градиентного спуска имеет свой параметр — скорость обучения. Обычно его подстаивают автоматически.
  8. Метод градиентного спуска повторяют много раз до тех пор, пока функция ошибки не перестанет значимо изменяться.

Регрессия с несколькими переменными

Множественная линейная регрессия

Множественная регрессия

Парная регрессия, как мы увидели выше, имеет дело с объектами, которые характеризуются одним числовым признаком ($x$). На практике, конечно, объекты характеризуются несколькими признаками, а значит в модели должна быть не одна входящая переменная, а несколько (или, что то же самое, вектор). Линейная регрессия с несколькими переменными также известна как «множественная линейная регрессия». Введем обозначения для уравнений, где мы можем иметь любое количество входных переменных:

$ x^{(i)} $- вектор-столбец всех значений признаков i-го обучающего примера;

$ x_j^{(i)} $ — значение j-го признака i-го обучающего примера;

$ x_j $ — вектор j-го признака всех обучающих примеров;

m — количество примеров в обучающей выборке;

n — количество признаков;

X — матрица признаков;

b — вектор параметров регрессии.

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

Для удобства примем, что $ x_0^{(i)} = 1 $ для всех $i$. Другими словами, мы ведем некий суррогатный признак, для всех объектов равный единице. Это никак не сказывается на самой функции гипотезы, это лишь условность обозначения, но это сильно упростит математические выкладки, особенно в матричной форме.

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

Общий вид модели множественной линейной регрессии:

[h_b(x) = b_0 + b_1 x_1 + b_2 x_2 + … + b_n x_n]

Или в матричной форме:

[h_b(x) = X cdot vec{b}]

Используя определение матричного умножения, наша многопараметрическая функция гипотезы может быть кратко представлена в виде: $h(x) = B X$.

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

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

Функция ошибки для множественной линейной регрессии:

[J(b) = frac{1}{2m} sum_{i=1}^{m} (h_b(x^{(i)}) — y^{(i)})^2]

Или в матричной форме:

[J(b) = frac{1}{2m} (X b — vec{y})^T (X b — vec{y})]

Обратите внимание, что мы специально не раскрываем выражение (h_b(x^{(i)})). Это нужно, чтобы подчеркнуть, что форма функции ошибки не зависит от функции гипотезы, она выражается через нее.

Теперь нам нужно взять производную этой функции ошибки. Здесь уже нужно знать производную самой функции гипотезы, так как:

[frac{partial}{partial b_i} J =
frac{1}{m} sum_{i=1}^{m} (h_b(x^{(i)}) — y^{(i)}) cdot frac{partial}{partial b_i} h_b(x^{(i)})]

В такой формулировке мы представляем частные производные функции ошибки (градиент) через частную производную функции гипотезы. Это так называемое моделенезависимое представление градиента. Ведь для этой формулы совершенно неважно, какой функцией будет наша гипотеза. Пока она является дифференцируемой, мы можем использовать градиент ее функции ошибки. Именно поэтому метод градиентного спуска работает с любыми аналитическими моделями, и нам не нужно каждый раз заново “переизобретать” математику градиентного спуска, адаптировать ее к каждой конкретной модели машинного обучения. Достаточно изучить этот метод один раз, в общей форме.

Метод градиентного спуска для множественной регрессии определяется следующими уравнениями:

повторять до сходимости:

[b_0 := b_0 — alpha frac{1}{m} sum_{i=1}^{m} (h_b(x^{(i)}) — y^{(i)}) cdot x_0^{(i)}]

[b_1 := b_1 — alpha frac{1}{m} sum_{i=1}^{m} (h_b(x^{(i)}) — y^{(i)}) cdot x_1^{(i)}]

[b_2 := b_2 — alpha frac{1}{m} sum_{i=1}^{m} (h_b(x^{(i)}) — y^{(i)}) cdot x_2^{(i)}]

[…]

Или в матричной форме:

[b := b — frac{alpha}{m} X^T (X b — vec{y})]

Выводы:

  1. Множественная регрессия очень похожа на парную, но с большим количеством признаков.
  2. Для удобства и однообразия, почти всегда обозначают $x_0 = 1$.
  3. Признаки образуют матрицу, поэтому уравнения множественной регрессии часто приводят в матричной форме, так короче.
  4. Алгоритм градиентного спуска для множественной регрессии точно такой же, как и для парной.

Нормализация признаков

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

В таком случае форма функции ошибки будет очень вытянутой. Это не проблема для математической формализации градиентного спуска — при достаточно малых $alpha$ метод все равно рано или поздно сходится. Проблема в практической реализации. Получается, что если выбрать скорость обучения выше определенного предела по самому компактному признаку, спуск разойдется. Значит, скорость обучения надо делать меньше. Но тогда в направлении второго признака спуск будет проходить слишком медленно. И получается, что градиентный спуск потребует гораздо больше итераций для завершения.

Эту проблему можно решить если изменить диапазоны входных данных, чтобы они выражались величинами примерно одного порядка. Это не позволит одному измерению численно доминировать над другим. На практике применяют несколько алгоритмов нормализации, самые распространенные из которых — минимаксная нормализация и стандартизация или z-оценки.

Минимаксная нормализация — это изменение входных данных по следующей формуле:

[x’ = frac{x — x_{min}}{x_{max} — x_{min}}]

После преобразования все значения будут лежать в диапазоне $x in [0; 1]$.

Z-оценки или стандартизация производится по формуле:

[x’ = frac{x — M[x]}{sigma_x}]

В таком случае данный признак приводится к стандартному распределению, то есть такому, у которого среднее 0, а дисперсия — 1.

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

Целевая переменная не нормируется.

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

Выводы:

  1. Нормализация нужна для ускорения метода градиентного спуска.
  2. Есть два основных метода нормализации — минимаксная и стандартизация.
  3. Параметры нормализации высчитываются по обучающей выборке.
  4. Нормализация встроена в большинство библиотечных методов.
  5. Некоторые методы более чувствительны к нормализации, чем другие.
  6. Нормализацию лучше сделать, чем не делать.

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

Нелинейная регрессия

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

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

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

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

Например, если наша функция гипотезы
$ hat{y} = h_b (x) = b_0 + b_1 x $,
то мы можем добавить еще один признак, основанный на $ x_1 $, получив квадратичную функцию

[hat{y} = h_b (x) = b_0 + b_1 x + b_2 x^2]

или кубическую функцию

[hat{y} = h_b (x) = b_0 + b_1 x + b_2 x^2 + b_3 x^3]

В кубической функции мы по сути ввели два новых признака:
$ x_2 = x^2, x_3 = x^3 $.
Точно таким же образом, мы можем создать, например, такую функцию:

[hat{y} = h_b (x) = b_0 + b_1 x + b_2 sqrt{x}]

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

И вот такое представление нелинейной функции как множественной линейной позволяет нам без изменений воспользоваться алгоритмом градиентного спуска для множественной линейной регрессии. Только вместо $ x_2, x_3, … , x_n $ нам нужно будет подставить соответствующие функции от $ x_1 $.

Полиномиальная регрессия
Источник: Wikimedia.

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

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

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

Для регрессии с двумя признаками.

Линейная модель (полином степени 1):

[h_b (x) = b_0 + b_1 x_1 + b_2 x_2]

Квадратичная модель (полином степени 2):

[h_b (x) = b_0 + b_1 x + b_2 x_2 + b_3 x_1^2 + b_4 x_2^2 + b_5 x_1 x_2]

Кубическая модель (полином степени 3):

[hat{y} = h_b (x) = b_0 + b_1 x_1 + b_2 x_2 + b_3 x_1^2 + b_4 x_2^2 + b_5 x_1 x_2 + b_6 x_1^3 + b_7 x_2^3 + b_7 x_1^2 x_2 + b_8 x_1 x_2^2]

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

Выводы:

  1. Данные в датасете не всегда располагаются так, что их хорошо может описывать линейная функция.
  2. Для описания нелинейных зависимостей нужна более сложная, нелинейная модель.
  3. Чтобы не изобретать алгоритм обучения заново, можно просто ввести в модель суррогатные признаки.
  4. Суррогатный признак — это новый признак, который считается из существующих атрибутов.
  5. Чаще всего используют полиномиальную регрессию — это когда в модель вводят полиномиальные признаки — степени существующих атрибутов.
  6. Обычно берут все комбинации факторов до какой-то определенной степени полинома.
  7. Полиномиальная регрессия может аппроксимировать любую функцию, нужно только подобрать степень полинома.
  8. Чем больше степень полиномиальной регрессии, тем она сложнее и универсальнее, но вычислительно сложнее (экспоненциально).

Практическое построение регрессии

В данном разделе мы посмотрим, как можно реализовать методы линейной регрессии на практике. Сначала мы попробуем создать алгоритм регрессии с нуля, а затем воспользуемся библиотечной функцией. Это поможет нам более полно понять, как работают модели машинного обучения в целом и в библиотеке sckikit-learn (самом популярном инструменте для создания и обучения моделей на языке программирования Python) в частности.

Для понимания данного раздела предполагаем, что читатель знаком с основами языка программирования Python. Нам понадобится знание его базового синтаксиса, немного — объектно-ориентированного программирования, немного — использования стандартных библиотек и модулей. Никаких продвинутых возможностей языка (типа метапрограммирования или декораторов) мы использовать не будем.

Как должны быть представлены данные для машинного обучения?

Применение любых моделей машинного обучения начинается с подготовки данных в необходимом формате. Для этого очень удобными для нас будут библиотеки numpy и pandas. Они практически всегда используются совместно с библиотекой sckikit-learn и другими инструментами машинного обучения. В первую очередь мы будем использовать numpy для создания массивов и операций с векторами и матрицами. Pandas нам понадобится для работы с табличными структурами — датасетами.

Если вы хотите самостоятельно задать в явном виде данные обучающей выборки, то нет ничего лучше использования обычных массивов ndarray. Обычно в одном массиве хранятся значения атрибутов — x, а в другом — значения целевой переменной — y.

1
2
3
4
5
6
7
8
9
10
11
import numpy as np

x = np.array([1.46, 1.13, -2.30, 1.74, 0.04, 
    -0.61, 0.32, -0.76, 0.58, -1.10, 
     0.87, 1.62, -0.53, -0.25, -1.07, 
    -0.38, -0.17, -0.32, -2.06, -0.88, ])

y = np.array([101.16, 78.44, -159.24, 120.72, 2.92, 
    -42.33, 22.07, -52.67, 40.32, -76.10, 
     59.88, 112.38, -36.54, -17.25, -74.24, 
    -26.57, -11.93, -22.31, -142.54, -60.74,])

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

1
2
3
4
5
x = np.array([
  [0, 1, 2, 3, 4],
  [5, 4, 9, 6, 3],
  [7.8, -0.1, 0.0, -2.14, 10.7],
  ])

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

Но чаще всего вы не будете задавать исходные данные явно. Практически всегда их приходится читать из каких-либо входных файлов. Удобнее всего это сделать при помощи библиотеки pandas вот так:

1
2
3
4
import pandas as pd

x = pd.read_csv('x.csv', index_col=0)
y = pd.read_csv('y.csv', index_col=0)

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

1
2
3
4
5
6
7
8
import pandas as pd

data = pd.read_csv('data.csv', index_col=0)

y = data.Y
y = data["Y"]

x = data.drop(["Y"])

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

Если вы используете pandas или numpy для формирования массивов данных, то получившиеся переменные будут разных типов — DataFrame или ndarray, соответственно. Но на дальнейшую работу это не повлияет, так как интерфейс работы с этими структурами данных очень похож. Например, неважно, какие именно массивы мы используем, их можно изобразить на графике вот так:

1
2
3
4
5
import maiplotlib.pyplot as plt

plt.figure()
plt.scatter(x, y)
plt.show()

Конечно, такая визуализация будет работать только в случае задачи парной регрессии. Если x многомерно, то простой график использовать не получится.

Давайте соберем весь наш код вместе:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
import numpy as np
import pandas as pd
import maiplotlib.pyplot as plt

# x = pd.read_csv('x.csv', index_col=0)
x = np.array([1.46, 1.13, -2.30, 1.74, 0.04, 
    -0.61, 0.32, -0.76, 0.58, -1.10, 
     0.87, 1.62, -0.53, -0.25, -1.07, 
    -0.38, -0.17, -0.32, -2.06, -0.88, ])

# y = pd.read_csv('y.csv', index_col=0)
y = np.array([101.16, 78.44, -159.24, 120.72, 2.92, 
    -42.33, 22.07, -52.67, 40.32, -76.10, 
     59.88, 112.38, -36.54, -17.25, -74.24, 
    -26.57, -11.93, -22.31, -142.54, -60.74,])

plt.figure()
plt.scatter(x, y)
plt.show()

Это код генерирует вот такой вот график:

Данные для регрессии

Как работает метод машинного обучения “на пальцах”?

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

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

1
2
3
4
5
class hypothesis(object):
    """Модель парной линейной регрессии"""
    def __init__(self):
        self.b0 = 0
        self.b1 = 0

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

Реализуем метод, который принимает значение входной переменной и возвращает теоретическое значение выходной — это прямое действие нашей регрессии — метод предсказания результата по факторам (в случае парной регрессии — по одному фактору):

1
2
    def predict(self, x):
        return self.b0 + self.b1 * x

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

Теперь зададим функцию ошибки:

1
2
    def error(self, X, Y):    
        return sum((self.predict(X) - Y)**2) / (2 * len(X)) 

В данном случае мы используем простую функцию ошибки — среднеквадратическое отклонение (mean squared error, MSE). Можно использовать и другие функции ошибки. Именно вид функции ошибки будет определять то, какой вид регрессии мы реализуем. Существует много разных вариаций простого алгоритма регрессии. О большинстве распространенных методах регрессии можно почитать в официальной документации sklearn.

Теперь реализуем метод градиентного спуска. Он должен принимать массив X и массив Y и обновлять параметры регрессии в соответствии в формулами градиентного спуска:

1
2
3
4
5
6
    def BGD(self, X, Y):  
        alpha = 0.5
        dJ0 = sum(self.predict(X) - Y) /len(X)
        dJ1 = sum((self.predict(X) - Y) * X) /len(X)
        self.b0 -= alpha * dJ0
        self.b1 -= alpha * dJ1

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

Давайте создадим объект регрессии и проверим начальное значение ошибки. В примерах приведены значения на модельном наборе данных, но этот метод можно использовать на любых данных, которые подходят по формату — x и y должны быть одномерными массивами чисел.

1
2
3
4
5
6
7
8
hyp = hypothesis()
print(hyp.predict(0))
print(hyp.predict(100))
J = hyp.error(x, y)
print("initial error:", J)
0 
0 
initial error: 36271.58344889084

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

Теперь все готово к запуску градиентного спуска.

1
2
3
4
5
6
7
8
9
10
hyp.BGD(x, y)
J = hyp.error(x, y)
print("error after gradient descent:", J)
error after gradient descent: 6734.135540194945
X0 = np.linspace(60, 180, 100)
Y0 = hyp.predict(X0)
plt.figure()
plt.scatter(x, y)
plt.plot(X0, Y0, 'r')
plt.show()

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

1
2
3
4
5
6
7
8
9
10
11
12
13
    def BGD(self, X, Y, alpha=0.5, accuracy=0.01, max_steps=5000):
        step = 0        
        old_err = hyp.error(X, Y)
        new_err = hyp.error(X, Y)
        dJ = 1
        while (dJ > accuracy) and (step < max_steps):
            dJ0 = sum(self.predict(X) - Y) /len(X)
            dJ1 = sum((self.predict(X) - Y) * X) /len(X)
            self.b0 -= alpha * dJ0
            self.b1 -= alpha * dJ1            
            old_err = new_err
            new_err = hyp.error(X, Y)
            dJ = abs(old_err - new_err) 

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

Запустим наш градиентный спуск:

1
2
3
4
5
hyp = hypothesis()
hyp.BGD(x, y)
J = hyp.error(x, y)
print("error after gradient descent:", J)
error after gradient descent: 298.76881676471504

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

Посмотрим, как теперь наша регрессия выглядит на графике:

1
2
3
4
5
6
X0 = np.linspace(60, 180, 100)
Y0 = hyp.predict(X0)
plt.figure()
plt.scatter(x, y)
plt.plot(X0, Y0, 'r')
plt.show()

Обученная регрессия

Уже значительно лучше. Линия регрессии довольно похожа на оптимальную. Так ли это на самом деле, глядя на график, сказать сложно, для этого нужно проанализировать, как ошибка регрессии менялась со временем:

Как оценить качество регрессионной модели?

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
    def BGD(self, X, Y, alpha=0.1, accuracy=0.01, max_steps=1000):
        steps, errors = [], []
        step = 0        
        old_err = hyp.error(X, Y)
        new_err = hyp.error(X, Y) - 1
        dJ = 1
        while (dJ > accuracy) and (step < max_steps):
            dJ0 = sum(self.predict(X) - Y) /len(X)
            dJ1 = sum((self.predict(X) - Y) * X) /len(X)
            self.b0 -= alpha * dJ0
            self.b1 -= alpha * dJ1            
            old_err = new_err
            new_err = hyp.error(X, Y)
            dJ = abs(old_err - new_err) 
            step += 1            
            steps.append(step)
            errors.append(new_err)
        return steps, errors

Мы просто запоминаем в массивах на номер шаа и ошибку на каждом шаге. Получив эти данные можно легко построить их на графике:

1
2
3
4
5
6
hyp = hypothesis()
steps, errors = hyp.BGD(x, y)

plt.figure()
plt.plot(steps, errors, 'g')
plt.show()

Прогресс обучения

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

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

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

Как подбирать скорость обучения?

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

На самом деле подобрать скорость обучения гораздо легче. Нужно использовать тот факт, что при превышении определенного порогового значения ошибка начинает возрастать. Кроме того, мы знаем, что скорость обучения должна быть положительна, но меньше единицы. Вся проблема в этом пороговом значении, которое сильно зависит от размерности задачи. При одних данных хорошо работает $ alpha = 0.5 $, а при каких-то приходится уменьшать ее на несколько порядков, например, $ alpha = 0.00000001 $.

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

Подбирать скорость обучения можно по следующему алгоритму. Сначала мы выбираем $ alpha $ близкое к 1, скажем, $ alpha = 0.7 $. Производим одну итерацию градиентного спуска и оцениваем, как изменилась ошибка. Если она уменьшилась, то ничего не надо менять, продолжаем спуск как обычно. Если же ошибка увеличилась, то скорость обучения нужно уменьшить. Например, раа в два. После чего мы повторяем первый шаг градиентного спуска. Таким образом мы не начинаем спуск, пока скорость обучения не снизится настолько, чтобы он начал сходиться.

Как применять регрессию с использованием scikit-learn?

Для серьезной работы, все-таки рекомендуется использовать готовые библиотечные решения. Они работаю гораздо быстрее, надежнее и гораздо проще, чем написанные самостоятельно. Мы будем использовать библиотеку scikit-learn для языка программирования Python как наш основной инструмент реализации простых моделей. Сегодня это одна их самых популярных библиотек для машинного обучения. Мы не будем повторять официальную документацию этой библиотеки, которая на редкость подробная и понятная. Наша задача — на примере этих инструментов понять, как работают и как применяются модели машинного обучения.

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

1
from sklearn import linear_model

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

Если вы используете DataFrame, то они обычно всегда настроены правильно, поэтому этого шага может не потребоваться. Важно запомнить, что все методы библиотечных моделей машинного обучения предполагают, что в x будет двумерный массив или DataFrame, а в y, соответственно, одномерный массив или Series.

Эта строка преобразует любой массив в вектор-столбец. Это если у вас один признак, то есть парная регрессия. Если признаков несколько, то вместо 1 следует указать число признаков. -1 на первой позиции означает, что по нулевому измерению будет столько элементов, сколько останется в массиве.

Само использование модели машинного обучения в этой библиотеке очень просто и сводится к трем действиям: создание экземпляра модели, обучение модели методом fit(), получение предсказаний методом predict(). Это общее поведение для любых моделей библиотеки. Для модели парной линейной регрессии нам понадобится класс LinearRegression.

1
2
3
4
5
6
reg = linear_model.LinearRegression()
reg.fit(x, y)
y_pred = reg.predict(x)

print(reg.score(x, y))
print("Коэффициенты: n", reg.coef_)

В этом классе кроме уже упомянутых методов fit() и predict(), которые есть в любой модели, есть большое количество методов и полей для получения дополнительной информации о моделях. Так, практически в каждой модели есть встроенный метод score(), который оценивает качество полученной модели. А поле coef_ содержит коэффициенты модели.

Обратите внимание, что в большинстве моделей коэффициентами считаются именно параметры при входящих переменных, то есть $ b_1, b_2, …, b_n $. Коэффициент $b_0$ считается особым и хранится отдельно в поле intercept_

Так как мы работаем с парной линейной регрессией, результат можно нарисовать на графике:

1
2
3
4
plt.figure(figsize=(12, 9))
plt.scatter(x, y, color="black")
plt.plot(x, y_pred, color="blue", linewidth=3)
plt.show()

Как мы видим, результат ничем не отличается от модели, которую мы обучили сами, вручную:

Библиотечная регрессия

Соберем код вместе и получим пример довольно реалистичного фрагмента работы с моделью машинного обучение. Примерно такой код можно встретить и в промышленных проектах по интеллектуальному анализу данных:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
from sklearn.linear_model import LinearRegression

x = x.reshape((-1, 1))

reg = LinearRegression()
reg.fit(x, y)
print(reg.score(x, y))

from sklearn.metrics import mean_squared_error, r2_score

y_pred = reg.predict(x)
print("Коэффициенты: n", reg.coef_)
print("Среднеквадратичная ошибка: %.2f" % mean_squared_error(y, y_pred))
print("Коэффициент детерминации: %.2f" % r2_score(y, y_pred))

plt.figure(figsize=(12, 9))
plt.scatter(x, y, color="black")
plt.plot(x, y_pred, color="blue", linewidth=3)
plt.show()


  Перевод


  Ссылка на автора

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

  • Тип задачи: регрессия? Классификация?
  • Бизнес цель?
  • Каково распределение целевой переменной?

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

Остров Андрос в Греции

Метрики регрессии

  • Средняя квадратическая ошибка (MSE)
  • Среднеквадратическая ошибка (RMSE)
  • Средняя абсолютная ошибка (MAE)
  • R в квадрате (R²)
  • Скорректированный R квадрат (R²)
  • Среднеквадратичная ошибка в процентах (MSPE)
  • Средняя абсолютная ошибка в процентах (MAPE)
  • Среднеквадратичная логарифмическая ошибка (RMSLE)

Средняя квадратическая ошибка (MSE)

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

гдеyᵢфактический ожидаемый результат иŷᵢэто прогноз модели.

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

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

Преимущество:Полезно, если у нас есть неожиданные значения, о которых мы должны заботиться. Очень высокое или низкое значение, на которое мы должны обратить внимание.

Недостаток:Если мы сделаем один очень плохой прогноз, возведение в квадрат сделает ошибку еще хуже, и это может исказить метрику в сторону переоценки плохости модели. Это особенно проблематичное поведение, если у нас есть зашумленные данные (то есть данные, которые по какой-либо причине не совсем надежны) — даже в «идеальной» модели может быть высокий MSE в этой ситуации, поэтому становится трудно судить, насколько хорошо модель выполняет. С другой стороны, если все ошибки малы или, скорее, меньше 1, то ощущается противоположный эффект: мы можем недооценивать недостатки модели.

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

Среднеквадратическая ошибка (RMSE)

RMSE — это просто квадратный корень из MSE. Квадратный корень введен, чтобы масштаб ошибок был таким же, как масштаб целей.

Теперь очень важно понять, в каком смысле RMSE похож на MSE, и в чем разница.

Во-первых, они похожи с точки зрения их минимизаторов, каждый минимизатор MSE также является минимизатором для RMSE и наоборот, поскольку квадратный корень является неубывающей функцией. Например, если у нас есть два набора предсказаний, A и B, и скажем, что MSE для A больше, чем MSE для B, то мы можем быть уверены, что RMSE для A больше RMSE для B. И это также работает в противоположном направлении. ,

Что это значит для нас?

Это означает, что, если целевым показателем является RMSE, мы все равно можем сравнивать наши модели, используя MSE, поскольку MSE упорядочит модели так же, как RMSE. Таким образом, мы можем оптимизировать MSE вместо RMSE.

На самом деле, с MSE работать немного проще, поэтому все используют MSE вместо RMSE. Также есть небольшая разница между этими двумя моделями на основе градиента.

Градиент СКО относительно i-го прогноза

Это означает, что путешествие по градиенту MSE эквивалентно путешествию по градиенту RMSE, но с другой скоростью потока, и скорость потока зависит от самой оценки MSE.

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

Средняя абсолютная ошибка (MAE)

В MAE ошибка рассчитывается как среднее абсолютных разностей между целевыми значениями и прогнозами. MAE — это линейная оценка, которая означает, чтовсе индивидуальные различия взвешены одинаковов среднем. Например, разница между 10 и 0 будет вдвое больше разницы между 5 и 0. Однако то же самое не верно для RMSE. Математически он рассчитывается по следующей формуле:

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

MAE широко используется в финансах, где ошибка в 10 долларов обычно в два раза хуже, чем ошибка в 5 долларов. С другой стороны, метрика MSE считает, что ошибка в 10 долларов в четыре раза хуже, чем ошибка в 5 долларов. MAE легче обосновать, чем RMSE.

Еще одна важная вещь в MAE — это его градиенты относительно прогнозов. Gradiend — это пошаговая функция, которая принимает -1, когда Y_hat меньше цели, и +1, когда она больше.

Теперь градиент не определен, когда предсказание является совершенным, потому что, когда Y_hat равен Y, мы не можем оценить градиент. Это не определено.

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

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

R в квадрате (R²)

А что если я скажу вам, что MSE для моих моделей предсказаний составляет 32? Должен ли я улучшить свою модель или она достаточно хороша? Или что, если мой MSE был 0,4? На самом деле, трудно понять, хороша наша модель или нет, посмотрев на абсолютные значения MSE или RMSE. Мы, вероятно, захотим измерить, как Во многом наша модель лучше, чем постоянная базовая линия.

Коэффициент детерминации, или R² (иногда читаемый как R-два), является еще одним показателем, который мы можем использовать для оценки модели, и он тесно связан с MSE, но имеет преимущество в том, чтобезмасштабное— не имеет значения, являются ли выходные значения очень большими или очень маленькими,R² всегда будет между -∞ и 1.

Когда R² отрицательно, это означает, что модель хуже, чем предсказание среднего значения.

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

гдеYс чертой означает среднее из наблюдаемогоyᵢ.

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

В заключение, R² — это соотношение между тем, насколько хороша наша модель, и тем, насколько хороша модель наивного среднего.

Распространенное заблуждение:Многие статьи в Интернете утверждают, что диапазон R² лежит между 0 и 1, что на самом деле не соответствует действительности. Максимальное значение R² равно 1, но минимальное может быть минус бесконечность.

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

MAE против MSE

Я заявил, что MAE более устойчив (менее чувствителен к выбросам), чем MSE, но это не значит, что всегда лучше использовать MAE. Следующие вопросы помогут вам решить:

Взять домой сообщение

В этой статье мы обсудили несколько важных метрик регрессии. Сначала мы обсудили среднеквадратичную ошибку и поняли, что наилучшей константой для нее является среднее целевое значение. Среднеквадратичная ошибка и R² очень похожи на MSE с точки зрения оптимизации. Затем мы обсудили среднюю абсолютную ошибку и когда люди предпочитают использовать MAE вместо MSE.

Спасибо за чтение, и я с нетерпением жду, чтобы услышать ваши вопросы 🙂Наслаждайтесь!

P.SСледите за моей следующей статьей, которая изучает другие более продвинутые метрики регрессии. Если вы хотите больше узнать о мире машинного обучения, вы также можете подписаться на меня в Instagram, напишите мне напрямую или найди меня на linkedin, Я хотел бы услышать от вас.

Ресурсы:
[1] https://dmitryulyanov.github.io/about

Рассмотрим двумерную
случайную величину (X,Y),
где X
и Y

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

,

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

Функцию
называют
«наилучшим приближением» Y
в смысле метода
наименьших квадратов, если математическое
ожидание М [Yg(X)]2
принимает наименьшее
возможное значение; функцию g(x)
называют
среднеквадратической
регрессией
Y
на X.

Теорема. Линейная
средняя квадратическая регрессия
Y
на X
имеет вид

,

где тх=M(X),
т
y=M(Y),
,
,r =
μ
ху
/(σxσy)
коэффициент корреляции
величин
X
и
Y.

Доказательство. Введем в
рассмотрение функцию двух независимых
аргументов
и:

.

Учитывая, что М
(X—тх)
(Ymy
)= 0, M[(X—mx)*(Ymy)]=
μ
ху=r
σxσy,

и выполнив выкладки, получим

.

Исследуем функцию на
экстремум, для чего приравняем нулю
частные производные:

Отсюда

,.

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

Итак, линейная средняя
квадратическая регрессия Y
и X
имеет вид

,

или

.

Коэффициент
называюткоэффициентом регрессии Y
на
X, а прямую

(**)

называют прямой
среднеквадратической регрессии
Y
на X.

Подставив найденные значения
ив соотношение (*), получим минимальное
значение функцииF
(,),
равное(1
r2).
Величину
(1
r2)
называют остаточной
дисперсией
случайной
величины Y
относительно случайной величины X;
она характеризует
величину ошибки, которую допускают при
замене Y
линейной функцией
g(X)=+X.
При r
= ±1 остаточная дисперсия равна нулю;
другими словами, при этих крайних
значениях коэффициента корреляции не
возникает ошибки при представлении Y
в виде линейной
функции от X.

Итак, если коэффициент
корреляции r
= ± 1, то Y и
X
связаны линейной
функциональной зависимостью.

Аналогично можно получить
прямую среднеквадратической регрессии
X
на Y:

(***)

(— коэффициент регрессии X
на Y)
и остаточную дисперсию(1
r2)величины
X
относительно Y.

Если r
= ±1, то обе прямые
регрессии, как видно из (**) и (***), совпадают.

Из уравнений (**) и (***) следует,
что обе прямые регрессии проходят через
точку (тх,
mу),
которую называют центром
совместного распределения величин
X
и Y.

§ 21. Линейная корреляция. Нормальная корреляция

Рассмотрим двумерную
случайную величину (X,
Y).
Если обе функции
регрессии Y
на X
и X
на Y
линейны (см. § 15), то говорят, что X
и Y связаны
линейной корреляционной
зависимостью.
Очевидно,
что графики линейных функций регрессии
— прямые линии, причем можно доказать,
что они совпадают с прямыми
среднеквадратической регрессии (см. §
20). Имеет место следующая важная теорема.

Теорема.
Если двумерная случайная
величина
(X,
Y)
распределена
нормально, то
X
и
Y связаны
линейной корреляционной зависимостью.

Доказательство. Двумерная плотность
вероятности (см. § 19)

,
(*)

где
,
.
(**)

Плотность вероятности
составляющей X
(см. § 19, замечание)

.
(***)

Найдем функцию регрессии
М (Y|х),
для чего сначала
найдем условный закон распределения
величины Y
при Х = х [см.
§ 14, формула (**)]:

.

Подставив (*) и (**) в правую часть этой
формулы выполнив выкладки, имеем

.

Заменив и
и v
по формулам (**),
окончательно получим

.

Полученное условное
распределение нормально с математическим
ожиданием (функцией регрессии Y
на X)

и дисперсией
.

Аналогично можно получить
функцию регрессии X
на Y:

.

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

Принимая во внимание вероятностный
смысл параметров двумерного нормального
распределения (см. § 19), заключаем, что
уравнения прямых регрессии

,

совпадают с
уравнениями прямых среднеквадратической
регрессии (см. § 20).

Задачи

1. Найти
законы распределения составляющих
дискретной двумерной случайной величины,
заданной законом распределения

Y

X

x1

x2

x3

y1

y2

0,12

0,10

0,18

0,11

0.10

0,39

Отв.
X x
1

x
2
x
3
Y y
1
y
2

p
0,22 0,29 0,49 p

0,40
0,60

2. Найти
вероятность того, что составляющая X
двумерной
случайной величины примет значение X
<
1/2
и
при этом составляющая Y
примет значение Y
<
1/3, если известна функция распределения
системы

.

Отв.
Р
(Х
<
1/2,
Y
<
1/3)
=
9/16.

3. Найти
вероятность попадания случайной
точки (X;
Y)
в
прямоугольник, ограниченный пряными
x=π/4,
x
= π/2, у = π/6, у = π/3,
если известна
функция распределения

F
(х,
у
)
=
sin
x
sin
у
(0
х
≤ π
/2,
0
y
π/2).

Отв.
P
(π/4
<
X
<
π/
2,
π/6
<
Y
<
π
/3)
=
0,11.

4. Найти плотность
распределения системы двух случайных
величин по известной функции распределения

F
(х,
у)=(12x)(13y)
(x0,
y0).

Отв..

5.
Внутри прямоугольника, ограниченного
прямыми x
=
0,
x
= π/2, у
=
0, y
=
π/2, плотность распределения системы
двух случайных величин f(x,
,y)=C
sin
(x
+ y
);
вне
прямоугольника f(x,
y)
=
0.
Найти:
а) величину С; б) функцию распределения
системы.

Отв.
а) С = 0,5; б) F(x,
у)
=0,5 [sin
x + sin у—sin
(x
+
y)]
(0xπ/2,
0yπ/2).

6.
Система двух случайных величин
распределена равномерно: в прямоугольнике,
ограниченном прямыми х
=
4,
х
= 6, y=10,
у
=
15, функция f
(x,
у)
сохраняет
постоянное значение, а вне этого
прямоугольника она равна нулю. Найти:
а) плотность f
(х,
у
)
совместного
распределения; б) функцию распределения
системы.

Отв.
а
)
f(x,y)=




б)
F(x,y)=.

7.
Плотность совместного распределения
системы двух случайных величин

f(x,y)=.

Найти:
а) величину C;
б) функцию распределения системы.

Отв.
а
)C=62
;
б
)
F(x,y)=.

8. Двумерная
случайная величина задана плотностью
совместного распределения

.

Найти
условные законы
распределения
составляющих.

Отв.
;.

ЧАСТЬ ТРЕТЬЯ

ЭЛЕМЕНТЫ МАТЕМАТИЧЕСКОЙ
СТАТИСТИКИ

Глава пятнадцатая

ВЫБОРОЧНЫЙ
МЕТОД

Соседние файлы в предмете [НЕСОРТИРОВАННОЕ]

  • #
  • #
  • #
  • #
  • #
  • #
  • #
  • #
  • #
  • #
  • #

In statistics, the mean squared error (MSE)[1] or mean squared deviation (MSD) of an estimator (of a procedure for estimating an unobserved quantity) measures the average of the squares of the errors—that is, the average squared difference between the estimated values and the actual value. MSE is a risk function, corresponding to the expected value of the squared error loss.[2] The fact that MSE is almost always strictly positive (and not zero) is because of randomness or because the estimator does not account for information that could produce a more accurate estimate.[3] In machine learning, specifically empirical risk minimization, MSE may refer to the empirical risk (the average loss on an observed data set), as an estimate of the true MSE (the true risk: the average loss on the actual population distribution).

The MSE is a measure of the quality of an estimator. As it is derived from the square of Euclidean distance, it is always a positive value that decreases as the error approaches zero.

The MSE is the second moment (about the origin) of the error, and thus incorporates both the variance of the estimator (how widely spread the estimates are from one data sample to another) and its bias (how far off the average estimated value is from the true value).[citation needed] For an unbiased estimator, the MSE is the variance of the estimator. Like the variance, MSE has the same units of measurement as the square of the quantity being estimated. In an analogy to standard deviation, taking the square root of MSE yields the root-mean-square error or root-mean-square deviation (RMSE or RMSD), which has the same units as the quantity being estimated; for an unbiased estimator, the RMSE is the square root of the variance, known as the standard error.

Definition and basic properties[edit]

The MSE either assesses the quality of a predictor (i.e., a function mapping arbitrary inputs to a sample of values of some random variable), or of an estimator (i.e., a mathematical function mapping a sample of data to an estimate of a parameter of the population from which the data is sampled). The definition of an MSE differs according to whether one is describing a predictor or an estimator.

Predictor[edit]

If a vector of n predictions is generated from a sample of n data points on all variables, and Y is the vector of observed values of the variable being predicted, with hat{Y} being the predicted values (e.g. as from a least-squares fit), then the within-sample MSE of the predictor is computed as

{displaystyle operatorname {MSE} ={frac {1}{n}}sum _{i=1}^{n}left(Y_{i}-{hat {Y_{i}}}right)^{2}.}

In other words, the MSE is the mean {textstyle left({frac {1}{n}}sum _{i=1}^{n}right)} of the squares of the errors {textstyle left(Y_{i}-{hat {Y_{i}}}right)^{2}}. This is an easily computable quantity for a particular sample (and hence is sample-dependent).

In matrix notation,

{displaystyle operatorname {MSE} ={frac {1}{n}}sum _{i=1}^{n}(e_{i})^{2}={frac {1}{n}}mathbf {e} ^{mathsf {T}}mathbf {e} }

where e_{i} is {displaystyle (Y_{i}-{hat {Y_{i}}})} and {displaystyle mathbf {e} } is the {displaystyle ntimes 1} column vector.

The MSE can also be computed on q data points that were not used in estimating the model, either because they were held back for this purpose, or because these data have been newly obtained. Within this process, known as statistical learning, the MSE is often called the test MSE,[4] and is computed as

{displaystyle operatorname {MSE} ={frac {1}{q}}sum _{i=n+1}^{n+q}left(Y_{i}-{hat {Y_{i}}}right)^{2}.}

Estimator[edit]

The MSE of an estimator hat{theta} with respect to an unknown parameter theta is defined as[1]

{displaystyle operatorname {MSE} ({hat {theta }})=operatorname {E} _{theta }left[({hat {theta }}-theta )^{2}right].}

This definition depends on the unknown parameter, but the MSE is a priori a property of an estimator. The MSE could be a function of unknown parameters, in which case any estimator of the MSE based on estimates of these parameters would be a function of the data (and thus a random variable). If the estimator hat{theta} is derived as a sample statistic and is used to estimate some population parameter, then the expectation is with respect to the sampling distribution of the sample statistic.

The MSE can be written as the sum of the variance of the estimator and the squared bias of the estimator, providing a useful way to calculate the MSE and implying that in the case of unbiased estimators, the MSE and variance are equivalent.[5]

{displaystyle operatorname {MSE} ({hat {theta }})=operatorname {Var} _{theta }({hat {theta }})+operatorname {Bias} ({hat {theta }},theta )^{2}.}

Proof of variance and bias relationship[edit]

{displaystyle {begin{aligned}operatorname {MSE} ({hat {theta }})&=operatorname {E} _{theta }left[({hat {theta }}-theta )^{2}right]\&=operatorname {E} _{theta }left[left({hat {theta }}-operatorname {E} _{theta }[{hat {theta }}]+operatorname {E} _{theta }[{hat {theta }}]-theta right)^{2}right]\&=operatorname {E} _{theta }left[left({hat {theta }}-operatorname {E} _{theta }[{hat {theta }}]right)^{2}+2left({hat {theta }}-operatorname {E} _{theta }[{hat {theta }}]right)left(operatorname {E} _{theta }[{hat {theta }}]-theta right)+left(operatorname {E} _{theta }[{hat {theta }}]-theta right)^{2}right]\&=operatorname {E} _{theta }left[left({hat {theta }}-operatorname {E} _{theta }[{hat {theta }}]right)^{2}right]+operatorname {E} _{theta }left[2left({hat {theta }}-operatorname {E} _{theta }[{hat {theta }}]right)left(operatorname {E} _{theta }[{hat {theta }}]-theta right)right]+operatorname {E} _{theta }left[left(operatorname {E} _{theta }[{hat {theta }}]-theta right)^{2}right]\&=operatorname {E} _{theta }left[left({hat {theta }}-operatorname {E} _{theta }[{hat {theta }}]right)^{2}right]+2left(operatorname {E} _{theta }[{hat {theta }}]-theta right)operatorname {E} _{theta }left[{hat {theta }}-operatorname {E} _{theta }[{hat {theta }}]right]+left(operatorname {E} _{theta }[{hat {theta }}]-theta right)^{2}&&operatorname {E} _{theta }[{hat {theta }}]-theta ={text{const.}}\&=operatorname {E} _{theta }left[left({hat {theta }}-operatorname {E} _{theta }[{hat {theta }}]right)^{2}right]+2left(operatorname {E} _{theta }[{hat {theta }}]-theta right)left(operatorname {E} _{theta }[{hat {theta }}]-operatorname {E} _{theta }[{hat {theta }}]right)+left(operatorname {E} _{theta }[{hat {theta }}]-theta right)^{2}&&operatorname {E} _{theta }[{hat {theta }}]={text{const.}}\&=operatorname {E} _{theta }left[left({hat {theta }}-operatorname {E} _{theta }[{hat {theta }}]right)^{2}right]+left(operatorname {E} _{theta }[{hat {theta }}]-theta right)^{2}\&=operatorname {Var} _{theta }({hat {theta }})+operatorname {Bias} _{theta }({hat {theta }},theta )^{2}end{aligned}}}

An even shorter proof can be achieved using the well-known formula that for a random variable {textstyle X}, {textstyle mathbb {E} (X^{2})=operatorname {Var} (X)+(mathbb {E} (X))^{2}}. By substituting {textstyle X} with, {textstyle {hat {theta }}-theta }, we have

{displaystyle {begin{aligned}operatorname {MSE} ({hat {theta }})&=mathbb {E} [({hat {theta }}-theta )^{2}]\&=operatorname {Var} ({hat {theta }}-theta )+(mathbb {E} [{hat {theta }}-theta ])^{2}\&=operatorname {Var} ({hat {theta }})+operatorname {Bias} ^{2}({hat {theta }})end{aligned}}}

But in real modeling case, MSE could be described as the addition of model variance, model bias, and irreducible uncertainty (see Bias–variance tradeoff). According to the relationship, the MSE of the estimators could be simply used for the efficiency comparison, which includes the information of estimator variance and bias. This is called MSE criterion.

In regression[edit]

In regression analysis, plotting is a more natural way to view the overall trend of the whole data. The mean of the distance from each point to the predicted regression model can be calculated, and shown as the mean squared error. The squaring is critical to reduce the complexity with negative signs. To minimize MSE, the model could be more accurate, which would mean the model is closer to actual data. One example of a linear regression using this method is the least squares method—which evaluates appropriateness of linear regression model to model bivariate dataset,[6] but whose limitation is related to known distribution of the data.

The term mean squared error is sometimes used to refer to the unbiased estimate of error variance: the residual sum of squares divided by the number of degrees of freedom. This definition for a known, computed quantity differs from the above definition for the computed MSE of a predictor, in that a different denominator is used. The denominator is the sample size reduced by the number of model parameters estimated from the same data, (np) for p regressors or (np−1) if an intercept is used (see errors and residuals in statistics for more details).[7] Although the MSE (as defined in this article) is not an unbiased estimator of the error variance, it is consistent, given the consistency of the predictor.

In regression analysis, «mean squared error», often referred to as mean squared prediction error or «out-of-sample mean squared error», can also refer to the mean value of the squared deviations of the predictions from the true values, over an out-of-sample test space, generated by a model estimated over a particular sample space. This also is a known, computed quantity, and it varies by sample and by out-of-sample test space.

Examples[edit]

Mean[edit]

Suppose we have a random sample of size n from a population, X_{1},dots ,X_{n}. Suppose the sample units were chosen with replacement. That is, the n units are selected one at a time, and previously selected units are still eligible for selection for all n draws. The usual estimator for the mu is the sample average

overline{X}=frac{1}{n}sum_{i=1}^n X_i

which has an expected value equal to the true mean mu (so it is unbiased) and a mean squared error of

{displaystyle operatorname {MSE} left({overline {X}}right)=operatorname {E} left[left({overline {X}}-mu right)^{2}right]=left({frac {sigma }{sqrt {n}}}right)^{2}={frac {sigma ^{2}}{n}}}

where sigma ^{2} is the population variance.

For a Gaussian distribution, this is the best unbiased estimator (i.e., one with the lowest MSE among all unbiased estimators), but not, say, for a uniform distribution.

Variance[edit]

The usual estimator for the variance is the corrected sample variance:

{displaystyle S_{n-1}^{2}={frac {1}{n-1}}sum _{i=1}^{n}left(X_{i}-{overline {X}}right)^{2}={frac {1}{n-1}}left(sum _{i=1}^{n}X_{i}^{2}-n{overline {X}}^{2}right).}

This is unbiased (its expected value is sigma ^{2}), hence also called the unbiased sample variance, and its MSE is[8]

{displaystyle operatorname {MSE} (S_{n-1}^{2})={frac {1}{n}}left(mu _{4}-{frac {n-3}{n-1}}sigma ^{4}right)={frac {1}{n}}left(gamma _{2}+{frac {2n}{n-1}}right)sigma ^{4},}

where mu _{4} is the fourth central moment of the distribution or population, and gamma_2=mu_4/sigma^4-3 is the excess kurtosis.

However, one can use other estimators for sigma ^{2} which are proportional to S^2_{n-1}, and an appropriate choice can always give a lower mean squared error. If we define

{displaystyle S_{a}^{2}={frac {n-1}{a}}S_{n-1}^{2}={frac {1}{a}}sum _{i=1}^{n}left(X_{i}-{overline {X}},right)^{2}}

then we calculate:

{displaystyle {begin{aligned}operatorname {MSE} (S_{a}^{2})&=operatorname {E} left[left({frac {n-1}{a}}S_{n-1}^{2}-sigma ^{2}right)^{2}right]\&=operatorname {E} left[{frac {(n-1)^{2}}{a^{2}}}S_{n-1}^{4}-2left({frac {n-1}{a}}S_{n-1}^{2}right)sigma ^{2}+sigma ^{4}right]\&={frac {(n-1)^{2}}{a^{2}}}operatorname {E} left[S_{n-1}^{4}right]-2left({frac {n-1}{a}}right)operatorname {E} left[S_{n-1}^{2}right]sigma ^{2}+sigma ^{4}\&={frac {(n-1)^{2}}{a^{2}}}operatorname {E} left[S_{n-1}^{4}right]-2left({frac {n-1}{a}}right)sigma ^{4}+sigma ^{4}&&operatorname {E} left[S_{n-1}^{2}right]=sigma ^{2}\&={frac {(n-1)^{2}}{a^{2}}}left({frac {gamma _{2}}{n}}+{frac {n+1}{n-1}}right)sigma ^{4}-2left({frac {n-1}{a}}right)sigma ^{4}+sigma ^{4}&&operatorname {E} left[S_{n-1}^{4}right]=operatorname {MSE} (S_{n-1}^{2})+sigma ^{4}\&={frac {n-1}{na^{2}}}left((n-1)gamma _{2}+n^{2}+nright)sigma ^{4}-2left({frac {n-1}{a}}right)sigma ^{4}+sigma ^{4}end{aligned}}}

This is minimized when

a=frac{(n-1)gamma_2+n^2+n}{n} = n+1+frac{n-1}{n}gamma_2.

For a Gaussian distribution, where gamma_2=0, this means that the MSE is minimized when dividing the sum by a=n+1. The minimum excess kurtosis is gamma_2=-2,[a] which is achieved by a Bernoulli distribution with p = 1/2 (a coin flip), and the MSE is minimized for {displaystyle a=n-1+{tfrac {2}{n}}.} Hence regardless of the kurtosis, we get a «better» estimate (in the sense of having a lower MSE) by scaling down the unbiased estimator a little bit; this is a simple example of a shrinkage estimator: one «shrinks» the estimator towards zero (scales down the unbiased estimator).

Further, while the corrected sample variance is the best unbiased estimator (minimum mean squared error among unbiased estimators) of variance for Gaussian distributions, if the distribution is not Gaussian, then even among unbiased estimators, the best unbiased estimator of the variance may not be S^2_{n-1}.

Gaussian distribution[edit]

The following table gives several estimators of the true parameters of the population, μ and σ2, for the Gaussian case.[9]

True value Estimator Mean squared error
{displaystyle theta =mu } hat{theta} = the unbiased estimator of the population mean, overline{X}=frac{1}{n}sum_{i=1}^n(X_i) operatorname{MSE}(overline{X})=operatorname{E}((overline{X}-mu)^2)=left(frac{sigma}{sqrt{n}}right)^2
{displaystyle theta =sigma ^{2}} hat{theta} = the unbiased estimator of the population variance, S^2_{n-1} = frac{1}{n-1}sum_{i=1}^nleft(X_i-overline{X},right)^2 operatorname{MSE}(S^2_{n-1})=operatorname{E}((S^2_{n-1}-sigma^2)^2)=frac{2}{n - 1}sigma^4
{displaystyle theta =sigma ^{2}} hat{theta} = the biased estimator of the population variance, S^2_{n} = frac{1}{n}sum_{i=1}^nleft(X_i-overline{X},right)^2 operatorname{MSE}(S^2_{n})=operatorname{E}((S^2_{n}-sigma^2)^2)=frac{2n - 1}{n^2}sigma^4
{displaystyle theta =sigma ^{2}} hat{theta} = the biased estimator of the population variance, S^2_{n+1} = frac{1}{n+1}sum_{i=1}^nleft(X_i-overline{X},right)^2 operatorname{MSE}(S^2_{n+1})=operatorname{E}((S^2_{n+1}-sigma^2)^2)=frac{2}{n + 1}sigma^4

Interpretation[edit]

An MSE of zero, meaning that the estimator hat{theta} predicts observations of the parameter theta with perfect accuracy, is ideal (but typically not possible).

Values of MSE may be used for comparative purposes. Two or more statistical models may be compared using their MSEs—as a measure of how well they explain a given set of observations: An unbiased estimator (estimated from a statistical model) with the smallest variance among all unbiased estimators is the best unbiased estimator or MVUE (Minimum-Variance Unbiased Estimator).

Both analysis of variance and linear regression techniques estimate the MSE as part of the analysis and use the estimated MSE to determine the statistical significance of the factors or predictors under study. The goal of experimental design is to construct experiments in such a way that when the observations are analyzed, the MSE is close to zero relative to the magnitude of at least one of the estimated treatment effects.

In one-way analysis of variance, MSE can be calculated by the division of the sum of squared errors and the degree of freedom. Also, the f-value is the ratio of the mean squared treatment and the MSE.

MSE is also used in several stepwise regression techniques as part of the determination as to how many predictors from a candidate set to include in a model for a given set of observations.

Applications[edit]

  • Minimizing MSE is a key criterion in selecting estimators: see minimum mean-square error. Among unbiased estimators, minimizing the MSE is equivalent to minimizing the variance, and the estimator that does this is the minimum variance unbiased estimator. However, a biased estimator may have lower MSE; see estimator bias.
  • In statistical modelling the MSE can represent the difference between the actual observations and the observation values predicted by the model. In this context, it is used to determine the extent to which the model fits the data as well as whether removing some explanatory variables is possible without significantly harming the model’s predictive ability.
  • In forecasting and prediction, the Brier score is a measure of forecast skill based on MSE.

Loss function[edit]

Squared error loss is one of the most widely used loss functions in statistics[citation needed], though its widespread use stems more from mathematical convenience than considerations of actual loss in applications. Carl Friedrich Gauss, who introduced the use of mean squared error, was aware of its arbitrariness and was in agreement with objections to it on these grounds.[3] The mathematical benefits of mean squared error are particularly evident in its use at analyzing the performance of linear regression, as it allows one to partition the variation in a dataset into variation explained by the model and variation explained by randomness.

Criticism[edit]

The use of mean squared error without question has been criticized by the decision theorist James Berger. Mean squared error is the negative of the expected value of one specific utility function, the quadratic utility function, which may not be the appropriate utility function to use under a given set of circumstances. There are, however, some scenarios where mean squared error can serve as a good approximation to a loss function occurring naturally in an application.[10]

Like variance, mean squared error has the disadvantage of heavily weighting outliers.[11] This is a result of the squaring of each term, which effectively weights large errors more heavily than small ones. This property, undesirable in many applications, has led researchers to use alternatives such as the mean absolute error, or those based on the median.

See also[edit]

  • Bias–variance tradeoff
  • Hodges’ estimator
  • James–Stein estimator
  • Mean percentage error
  • Mean square quantization error
  • Mean square weighted deviation
  • Mean squared displacement
  • Mean squared prediction error
  • Minimum mean square error
  • Minimum mean squared error estimator
  • Overfitting
  • Peak signal-to-noise ratio

Notes[edit]

  1. ^ This can be proved by Jensen’s inequality as follows. The fourth central moment is an upper bound for the square of variance, so that the least value for their ratio is one, therefore, the least value for the excess kurtosis is −2, achieved, for instance, by a Bernoulli with p=1/2.

References[edit]

  1. ^ a b «Mean Squared Error (MSE)». www.probabilitycourse.com. Retrieved 2020-09-12.
  2. ^ Bickel, Peter J.; Doksum, Kjell A. (2015). Mathematical Statistics: Basic Ideas and Selected Topics. Vol. I (Second ed.). p. 20. If we use quadratic loss, our risk function is called the mean squared error (MSE) …
  3. ^ a b Lehmann, E. L.; Casella, George (1998). Theory of Point Estimation (2nd ed.). New York: Springer. ISBN 978-0-387-98502-2. MR 1639875.
  4. ^ Gareth, James; Witten, Daniela; Hastie, Trevor; Tibshirani, Rob (2021). An Introduction to Statistical Learning: with Applications in R. Springer. ISBN 978-1071614174.
  5. ^ Wackerly, Dennis; Mendenhall, William; Scheaffer, Richard L. (2008). Mathematical Statistics with Applications (7 ed.). Belmont, CA, USA: Thomson Higher Education. ISBN 978-0-495-38508-0.
  6. ^ A modern introduction to probability and statistics : understanding why and how. Dekking, Michel, 1946-. London: Springer. 2005. ISBN 978-1-85233-896-1. OCLC 262680588.{{cite book}}: CS1 maint: others (link)
  7. ^ Steel, R.G.D, and Torrie, J. H., Principles and Procedures of Statistics with Special Reference to the Biological Sciences., McGraw Hill, 1960, page 288.
  8. ^ Mood, A.; Graybill, F.; Boes, D. (1974). Introduction to the Theory of Statistics (3rd ed.). McGraw-Hill. p. 229.
  9. ^ DeGroot, Morris H. (1980). Probability and Statistics (2nd ed.). Addison-Wesley.
  10. ^ Berger, James O. (1985). «2.4.2 Certain Standard Loss Functions». Statistical Decision Theory and Bayesian Analysis (2nd ed.). New York: Springer-Verlag. p. 60. ISBN 978-0-387-96098-2. MR 0804611.
  11. ^ Bermejo, Sergio; Cabestany, Joan (2001). «Oriented principal component analysis for large margin classifiers». Neural Networks. 14 (10): 1447–1461. doi:10.1016/S0893-6080(01)00106-X. PMID 11771723.

In statistics, the mean squared error (MSE)[1] or mean squared deviation (MSD) of an estimator (of a procedure for estimating an unobserved quantity) measures the average of the squares of the errors—that is, the average squared difference between the estimated values and the actual value. MSE is a risk function, corresponding to the expected value of the squared error loss.[2] The fact that MSE is almost always strictly positive (and not zero) is because of randomness or because the estimator does not account for information that could produce a more accurate estimate.[3] In machine learning, specifically empirical risk minimization, MSE may refer to the empirical risk (the average loss on an observed data set), as an estimate of the true MSE (the true risk: the average loss on the actual population distribution).

The MSE is a measure of the quality of an estimator. As it is derived from the square of Euclidean distance, it is always a positive value that decreases as the error approaches zero.

The MSE is the second moment (about the origin) of the error, and thus incorporates both the variance of the estimator (how widely spread the estimates are from one data sample to another) and its bias (how far off the average estimated value is from the true value).[citation needed] For an unbiased estimator, the MSE is the variance of the estimator. Like the variance, MSE has the same units of measurement as the square of the quantity being estimated. In an analogy to standard deviation, taking the square root of MSE yields the root-mean-square error or root-mean-square deviation (RMSE or RMSD), which has the same units as the quantity being estimated; for an unbiased estimator, the RMSE is the square root of the variance, known as the standard error.

Definition and basic properties[edit]

The MSE either assesses the quality of a predictor (i.e., a function mapping arbitrary inputs to a sample of values of some random variable), or of an estimator (i.e., a mathematical function mapping a sample of data to an estimate of a parameter of the population from which the data is sampled). The definition of an MSE differs according to whether one is describing a predictor or an estimator.

Predictor[edit]

If a vector of n predictions is generated from a sample of n data points on all variables, and Y is the vector of observed values of the variable being predicted, with hat{Y} being the predicted values (e.g. as from a least-squares fit), then the within-sample MSE of the predictor is computed as

{displaystyle operatorname {MSE} ={frac {1}{n}}sum _{i=1}^{n}left(Y_{i}-{hat {Y_{i}}}right)^{2}.}

In other words, the MSE is the mean {textstyle left({frac {1}{n}}sum _{i=1}^{n}right)} of the squares of the errors {textstyle left(Y_{i}-{hat {Y_{i}}}right)^{2}}. This is an easily computable quantity for a particular sample (and hence is sample-dependent).

In matrix notation,

{displaystyle operatorname {MSE} ={frac {1}{n}}sum _{i=1}^{n}(e_{i})^{2}={frac {1}{n}}mathbf {e} ^{mathsf {T}}mathbf {e} }

where e_{i} is {displaystyle (Y_{i}-{hat {Y_{i}}})} and {displaystyle mathbf {e} } is the {displaystyle ntimes 1} column vector.

The MSE can also be computed on q data points that were not used in estimating the model, either because they were held back for this purpose, or because these data have been newly obtained. Within this process, known as statistical learning, the MSE is often called the test MSE,[4] and is computed as

{displaystyle operatorname {MSE} ={frac {1}{q}}sum _{i=n+1}^{n+q}left(Y_{i}-{hat {Y_{i}}}right)^{2}.}

Estimator[edit]

The MSE of an estimator hat{theta} with respect to an unknown parameter theta is defined as[1]

{displaystyle operatorname {MSE} ({hat {theta }})=operatorname {E} _{theta }left[({hat {theta }}-theta )^{2}right].}

This definition depends on the unknown parameter, but the MSE is a priori a property of an estimator. The MSE could be a function of unknown parameters, in which case any estimator of the MSE based on estimates of these parameters would be a function of the data (and thus a random variable). If the estimator hat{theta} is derived as a sample statistic and is used to estimate some population parameter, then the expectation is with respect to the sampling distribution of the sample statistic.

The MSE can be written as the sum of the variance of the estimator and the squared bias of the estimator, providing a useful way to calculate the MSE and implying that in the case of unbiased estimators, the MSE and variance are equivalent.[5]

{displaystyle operatorname {MSE} ({hat {theta }})=operatorname {Var} _{theta }({hat {theta }})+operatorname {Bias} ({hat {theta }},theta )^{2}.}

Proof of variance and bias relationship[edit]

{displaystyle {begin{aligned}operatorname {MSE} ({hat {theta }})&=operatorname {E} _{theta }left[({hat {theta }}-theta )^{2}right]\&=operatorname {E} _{theta }left[left({hat {theta }}-operatorname {E} _{theta }[{hat {theta }}]+operatorname {E} _{theta }[{hat {theta }}]-theta right)^{2}right]\&=operatorname {E} _{theta }left[left({hat {theta }}-operatorname {E} _{theta }[{hat {theta }}]right)^{2}+2left({hat {theta }}-operatorname {E} _{theta }[{hat {theta }}]right)left(operatorname {E} _{theta }[{hat {theta }}]-theta right)+left(operatorname {E} _{theta }[{hat {theta }}]-theta right)^{2}right]\&=operatorname {E} _{theta }left[left({hat {theta }}-operatorname {E} _{theta }[{hat {theta }}]right)^{2}right]+operatorname {E} _{theta }left[2left({hat {theta }}-operatorname {E} _{theta }[{hat {theta }}]right)left(operatorname {E} _{theta }[{hat {theta }}]-theta right)right]+operatorname {E} _{theta }left[left(operatorname {E} _{theta }[{hat {theta }}]-theta right)^{2}right]\&=operatorname {E} _{theta }left[left({hat {theta }}-operatorname {E} _{theta }[{hat {theta }}]right)^{2}right]+2left(operatorname {E} _{theta }[{hat {theta }}]-theta right)operatorname {E} _{theta }left[{hat {theta }}-operatorname {E} _{theta }[{hat {theta }}]right]+left(operatorname {E} _{theta }[{hat {theta }}]-theta right)^{2}&&operatorname {E} _{theta }[{hat {theta }}]-theta ={text{const.}}\&=operatorname {E} _{theta }left[left({hat {theta }}-operatorname {E} _{theta }[{hat {theta }}]right)^{2}right]+2left(operatorname {E} _{theta }[{hat {theta }}]-theta right)left(operatorname {E} _{theta }[{hat {theta }}]-operatorname {E} _{theta }[{hat {theta }}]right)+left(operatorname {E} _{theta }[{hat {theta }}]-theta right)^{2}&&operatorname {E} _{theta }[{hat {theta }}]={text{const.}}\&=operatorname {E} _{theta }left[left({hat {theta }}-operatorname {E} _{theta }[{hat {theta }}]right)^{2}right]+left(operatorname {E} _{theta }[{hat {theta }}]-theta right)^{2}\&=operatorname {Var} _{theta }({hat {theta }})+operatorname {Bias} _{theta }({hat {theta }},theta )^{2}end{aligned}}}

An even shorter proof can be achieved using the well-known formula that for a random variable {textstyle X}, {textstyle mathbb {E} (X^{2})=operatorname {Var} (X)+(mathbb {E} (X))^{2}}. By substituting {textstyle X} with, {textstyle {hat {theta }}-theta }, we have

{displaystyle {begin{aligned}operatorname {MSE} ({hat {theta }})&=mathbb {E} [({hat {theta }}-theta )^{2}]\&=operatorname {Var} ({hat {theta }}-theta )+(mathbb {E} [{hat {theta }}-theta ])^{2}\&=operatorname {Var} ({hat {theta }})+operatorname {Bias} ^{2}({hat {theta }})end{aligned}}}

But in real modeling case, MSE could be described as the addition of model variance, model bias, and irreducible uncertainty (see Bias–variance tradeoff). According to the relationship, the MSE of the estimators could be simply used for the efficiency comparison, which includes the information of estimator variance and bias. This is called MSE criterion.

In regression[edit]

In regression analysis, plotting is a more natural way to view the overall trend of the whole data. The mean of the distance from each point to the predicted regression model can be calculated, and shown as the mean squared error. The squaring is critical to reduce the complexity with negative signs. To minimize MSE, the model could be more accurate, which would mean the model is closer to actual data. One example of a linear regression using this method is the least squares method—which evaluates appropriateness of linear regression model to model bivariate dataset,[6] but whose limitation is related to known distribution of the data.

The term mean squared error is sometimes used to refer to the unbiased estimate of error variance: the residual sum of squares divided by the number of degrees of freedom. This definition for a known, computed quantity differs from the above definition for the computed MSE of a predictor, in that a different denominator is used. The denominator is the sample size reduced by the number of model parameters estimated from the same data, (np) for p regressors or (np−1) if an intercept is used (see errors and residuals in statistics for more details).[7] Although the MSE (as defined in this article) is not an unbiased estimator of the error variance, it is consistent, given the consistency of the predictor.

In regression analysis, «mean squared error», often referred to as mean squared prediction error or «out-of-sample mean squared error», can also refer to the mean value of the squared deviations of the predictions from the true values, over an out-of-sample test space, generated by a model estimated over a particular sample space. This also is a known, computed quantity, and it varies by sample and by out-of-sample test space.

Examples[edit]

Mean[edit]

Suppose we have a random sample of size n from a population, X_{1},dots ,X_{n}. Suppose the sample units were chosen with replacement. That is, the n units are selected one at a time, and previously selected units are still eligible for selection for all n draws. The usual estimator for the mu is the sample average

overline{X}=frac{1}{n}sum_{i=1}^n X_i

which has an expected value equal to the true mean mu (so it is unbiased) and a mean squared error of

{displaystyle operatorname {MSE} left({overline {X}}right)=operatorname {E} left[left({overline {X}}-mu right)^{2}right]=left({frac {sigma }{sqrt {n}}}right)^{2}={frac {sigma ^{2}}{n}}}

where sigma ^{2} is the population variance.

For a Gaussian distribution, this is the best unbiased estimator (i.e., one with the lowest MSE among all unbiased estimators), but not, say, for a uniform distribution.

Variance[edit]

The usual estimator for the variance is the corrected sample variance:

{displaystyle S_{n-1}^{2}={frac {1}{n-1}}sum _{i=1}^{n}left(X_{i}-{overline {X}}right)^{2}={frac {1}{n-1}}left(sum _{i=1}^{n}X_{i}^{2}-n{overline {X}}^{2}right).}

This is unbiased (its expected value is sigma ^{2}), hence also called the unbiased sample variance, and its MSE is[8]

{displaystyle operatorname {MSE} (S_{n-1}^{2})={frac {1}{n}}left(mu _{4}-{frac {n-3}{n-1}}sigma ^{4}right)={frac {1}{n}}left(gamma _{2}+{frac {2n}{n-1}}right)sigma ^{4},}

where mu _{4} is the fourth central moment of the distribution or population, and gamma_2=mu_4/sigma^4-3 is the excess kurtosis.

However, one can use other estimators for sigma ^{2} which are proportional to S^2_{n-1}, and an appropriate choice can always give a lower mean squared error. If we define

{displaystyle S_{a}^{2}={frac {n-1}{a}}S_{n-1}^{2}={frac {1}{a}}sum _{i=1}^{n}left(X_{i}-{overline {X}},right)^{2}}

then we calculate:

{displaystyle {begin{aligned}operatorname {MSE} (S_{a}^{2})&=operatorname {E} left[left({frac {n-1}{a}}S_{n-1}^{2}-sigma ^{2}right)^{2}right]\&=operatorname {E} left[{frac {(n-1)^{2}}{a^{2}}}S_{n-1}^{4}-2left({frac {n-1}{a}}S_{n-1}^{2}right)sigma ^{2}+sigma ^{4}right]\&={frac {(n-1)^{2}}{a^{2}}}operatorname {E} left[S_{n-1}^{4}right]-2left({frac {n-1}{a}}right)operatorname {E} left[S_{n-1}^{2}right]sigma ^{2}+sigma ^{4}\&={frac {(n-1)^{2}}{a^{2}}}operatorname {E} left[S_{n-1}^{4}right]-2left({frac {n-1}{a}}right)sigma ^{4}+sigma ^{4}&&operatorname {E} left[S_{n-1}^{2}right]=sigma ^{2}\&={frac {(n-1)^{2}}{a^{2}}}left({frac {gamma _{2}}{n}}+{frac {n+1}{n-1}}right)sigma ^{4}-2left({frac {n-1}{a}}right)sigma ^{4}+sigma ^{4}&&operatorname {E} left[S_{n-1}^{4}right]=operatorname {MSE} (S_{n-1}^{2})+sigma ^{4}\&={frac {n-1}{na^{2}}}left((n-1)gamma _{2}+n^{2}+nright)sigma ^{4}-2left({frac {n-1}{a}}right)sigma ^{4}+sigma ^{4}end{aligned}}}

This is minimized when

a=frac{(n-1)gamma_2+n^2+n}{n} = n+1+frac{n-1}{n}gamma_2.

For a Gaussian distribution, where gamma_2=0, this means that the MSE is minimized when dividing the sum by a=n+1. The minimum excess kurtosis is gamma_2=-2,[a] which is achieved by a Bernoulli distribution with p = 1/2 (a coin flip), and the MSE is minimized for {displaystyle a=n-1+{tfrac {2}{n}}.} Hence regardless of the kurtosis, we get a «better» estimate (in the sense of having a lower MSE) by scaling down the unbiased estimator a little bit; this is a simple example of a shrinkage estimator: one «shrinks» the estimator towards zero (scales down the unbiased estimator).

Further, while the corrected sample variance is the best unbiased estimator (minimum mean squared error among unbiased estimators) of variance for Gaussian distributions, if the distribution is not Gaussian, then even among unbiased estimators, the best unbiased estimator of the variance may not be S^2_{n-1}.

Gaussian distribution[edit]

The following table gives several estimators of the true parameters of the population, μ and σ2, for the Gaussian case.[9]

True value Estimator Mean squared error
{displaystyle theta =mu } hat{theta} = the unbiased estimator of the population mean, overline{X}=frac{1}{n}sum_{i=1}^n(X_i) operatorname{MSE}(overline{X})=operatorname{E}((overline{X}-mu)^2)=left(frac{sigma}{sqrt{n}}right)^2
{displaystyle theta =sigma ^{2}} hat{theta} = the unbiased estimator of the population variance, S^2_{n-1} = frac{1}{n-1}sum_{i=1}^nleft(X_i-overline{X},right)^2 operatorname{MSE}(S^2_{n-1})=operatorname{E}((S^2_{n-1}-sigma^2)^2)=frac{2}{n - 1}sigma^4
{displaystyle theta =sigma ^{2}} hat{theta} = the biased estimator of the population variance, S^2_{n} = frac{1}{n}sum_{i=1}^nleft(X_i-overline{X},right)^2 operatorname{MSE}(S^2_{n})=operatorname{E}((S^2_{n}-sigma^2)^2)=frac{2n - 1}{n^2}sigma^4
{displaystyle theta =sigma ^{2}} hat{theta} = the biased estimator of the population variance, S^2_{n+1} = frac{1}{n+1}sum_{i=1}^nleft(X_i-overline{X},right)^2 operatorname{MSE}(S^2_{n+1})=operatorname{E}((S^2_{n+1}-sigma^2)^2)=frac{2}{n + 1}sigma^4

Interpretation[edit]

An MSE of zero, meaning that the estimator hat{theta} predicts observations of the parameter theta with perfect accuracy, is ideal (but typically not possible).

Values of MSE may be used for comparative purposes. Two or more statistical models may be compared using their MSEs—as a measure of how well they explain a given set of observations: An unbiased estimator (estimated from a statistical model) with the smallest variance among all unbiased estimators is the best unbiased estimator or MVUE (Minimum-Variance Unbiased Estimator).

Both analysis of variance and linear regression techniques estimate the MSE as part of the analysis and use the estimated MSE to determine the statistical significance of the factors or predictors under study. The goal of experimental design is to construct experiments in such a way that when the observations are analyzed, the MSE is close to zero relative to the magnitude of at least one of the estimated treatment effects.

In one-way analysis of variance, MSE can be calculated by the division of the sum of squared errors and the degree of freedom. Also, the f-value is the ratio of the mean squared treatment and the MSE.

MSE is also used in several stepwise regression techniques as part of the determination as to how many predictors from a candidate set to include in a model for a given set of observations.

Applications[edit]

  • Minimizing MSE is a key criterion in selecting estimators: see minimum mean-square error. Among unbiased estimators, minimizing the MSE is equivalent to minimizing the variance, and the estimator that does this is the minimum variance unbiased estimator. However, a biased estimator may have lower MSE; see estimator bias.
  • In statistical modelling the MSE can represent the difference between the actual observations and the observation values predicted by the model. In this context, it is used to determine the extent to which the model fits the data as well as whether removing some explanatory variables is possible without significantly harming the model’s predictive ability.
  • In forecasting and prediction, the Brier score is a measure of forecast skill based on MSE.

Loss function[edit]

Squared error loss is one of the most widely used loss functions in statistics[citation needed], though its widespread use stems more from mathematical convenience than considerations of actual loss in applications. Carl Friedrich Gauss, who introduced the use of mean squared error, was aware of its arbitrariness and was in agreement with objections to it on these grounds.[3] The mathematical benefits of mean squared error are particularly evident in its use at analyzing the performance of linear regression, as it allows one to partition the variation in a dataset into variation explained by the model and variation explained by randomness.

Criticism[edit]

The use of mean squared error without question has been criticized by the decision theorist James Berger. Mean squared error is the negative of the expected value of one specific utility function, the quadratic utility function, which may not be the appropriate utility function to use under a given set of circumstances. There are, however, some scenarios where mean squared error can serve as a good approximation to a loss function occurring naturally in an application.[10]

Like variance, mean squared error has the disadvantage of heavily weighting outliers.[11] This is a result of the squaring of each term, which effectively weights large errors more heavily than small ones. This property, undesirable in many applications, has led researchers to use alternatives such as the mean absolute error, or those based on the median.

See also[edit]

  • Bias–variance tradeoff
  • Hodges’ estimator
  • James–Stein estimator
  • Mean percentage error
  • Mean square quantization error
  • Mean square weighted deviation
  • Mean squared displacement
  • Mean squared prediction error
  • Minimum mean square error
  • Minimum mean squared error estimator
  • Overfitting
  • Peak signal-to-noise ratio

Notes[edit]

  1. ^ This can be proved by Jensen’s inequality as follows. The fourth central moment is an upper bound for the square of variance, so that the least value for their ratio is one, therefore, the least value for the excess kurtosis is −2, achieved, for instance, by a Bernoulli with p=1/2.

References[edit]

  1. ^ a b «Mean Squared Error (MSE)». www.probabilitycourse.com. Retrieved 2020-09-12.
  2. ^ Bickel, Peter J.; Doksum, Kjell A. (2015). Mathematical Statistics: Basic Ideas and Selected Topics. Vol. I (Second ed.). p. 20. If we use quadratic loss, our risk function is called the mean squared error (MSE) …
  3. ^ a b Lehmann, E. L.; Casella, George (1998). Theory of Point Estimation (2nd ed.). New York: Springer. ISBN 978-0-387-98502-2. MR 1639875.
  4. ^ Gareth, James; Witten, Daniela; Hastie, Trevor; Tibshirani, Rob (2021). An Introduction to Statistical Learning: with Applications in R. Springer. ISBN 978-1071614174.
  5. ^ Wackerly, Dennis; Mendenhall, William; Scheaffer, Richard L. (2008). Mathematical Statistics with Applications (7 ed.). Belmont, CA, USA: Thomson Higher Education. ISBN 978-0-495-38508-0.
  6. ^ A modern introduction to probability and statistics : understanding why and how. Dekking, Michel, 1946-. London: Springer. 2005. ISBN 978-1-85233-896-1. OCLC 262680588.{{cite book}}: CS1 maint: others (link)
  7. ^ Steel, R.G.D, and Torrie, J. H., Principles and Procedures of Statistics with Special Reference to the Biological Sciences., McGraw Hill, 1960, page 288.
  8. ^ Mood, A.; Graybill, F.; Boes, D. (1974). Introduction to the Theory of Statistics (3rd ed.). McGraw-Hill. p. 229.
  9. ^ DeGroot, Morris H. (1980). Probability and Statistics (2nd ed.). Addison-Wesley.
  10. ^ Berger, James O. (1985). «2.4.2 Certain Standard Loss Functions». Statistical Decision Theory and Bayesian Analysis (2nd ed.). New York: Springer-Verlag. p. 60. ISBN 978-0-387-96098-2. MR 0804611.
  11. ^ Bermejo, Sergio; Cabestany, Joan (2001). «Oriented principal component analysis for large margin classifiers». Neural Networks. 14 (10): 1447–1461. doi:10.1016/S0893-6080(01)00106-X. PMID 11771723.

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

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

Понятия ковариации и корреляции

Для примера будем использовать набор данных Animals из пакета MASS, описание которого можно получить с помощью команды ?Animals:

> head(Animals)
                    body brain
Mountain beaver     1.35   8.1
Cow               465.00 423.0
Grey wolf          36.33 119.5
Goat               27.66 115.0
Guinea pig          1.04   5.5
Dipliodocus     11700.00  50.0

> nrow(Animals)
[1] 28

У нас имеется n=28 наблюдений и мы хотим проверить «силу взаимосвязи» (и есть ли такая взаимосвязь вообще) между весом мозга (переменная brain) и весом тела (переменная body). Хорошим правилом считается посмотреть на данные, с которыми мы работаем, поэтому отобразим точки на диаграмме рассеяния:

plot_animals <- function(animals, xlab='', ylab='', main_label='') {
    plot(animals$body, animals$brain,
         xlab=xlab, ylab=ylab, main=main_label,
         col='royal blue', lwd=2)
}

plot_animals(Animals, xlab='body weight', ylab='brain weight',
             main='Brain and Body Weights for 28 Species')

Из полученного графика может показаться, что между наблюдениями (парами x_i и y_i) нет никакой взаимосвязи, поэтому попробуем преобразовать наши данные к log-log шкале:

animals.log <- data.frame(body=log(Animals$body), brain=log(Animals$brain))

plot_animals(animals.log, xlab='log(body weight)', ylab='log(brain weight)',
             main='Brain and Body Weights for 28 Species')

Note

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

Теперь стало более очевидно, что имеется линейная взаимосвязь между наблюдениями. Возникает вопрос: «На сколько эта взаимосвязь сильная?» Чтобы ответить на этот вопрос мы можем посчитать так называемый коэффициент корреляции. Для этого сначала найдем выборочные средние и отобразим их соответствующими линиями на диаграмме рассеяния:

bar{y} = frac{sum_{i=1}^{n}y_{i}}{n}

bar{x} = frac{sum_{i=1}^{n}x_{i}}{n}

plot_animals(animals.log, xlab='log(body weight)', ylab='log(brain weight)',
             main = 'Brain and Body Weights for 28 Species')

# Номера четвертей
text(x=10, y=8, labels="1", col="black", cex=2)
text(x=10, y=0, labels="4", col="black", cex=2)
text(x=-3, y=8, labels="2", col="black", cex=2)
text(x=-3, y=0, labels="3", col="black", cex=2)

# Знаки в четвертях
text(x=7, y=6, labels="+/+", col="dark red", cex=2)
text(x=0, y=6, labels="+/-", col="dark red", cex=2)
text(x=7, y=2, labels="-/+", col="dark red", cex=2)
text(x=0, y=2, labels="-/-", col="dark red", cex=2)

(body.mean <- mean(animals.log$body))
[1] 3.771306
(brain.mean <- mean(animals.log$brain))
[1] 4.425446

abline(v=body.mean, col="dark red", cex=2)
abline(h=brain.mean, col="dark red", cex=2)

Две линии делят диаграмму на четыре квадранта. Для каждой i-ой точки на графике вычислим следующие величины:

  • y_{i} — bar{y} отклонение каждого наблюдения y_{i} от выборочного среднего;
  • x_{i} — bar{x} отклонение каждого наблюдения x_{i} от выборочного среднего;
  • произведение двух величин (y_{i} — bar{y})(x_{i} — bar{x}).
y_deviation <- animals.log$brain - brain.mean
x_deviation <- animals.log$body - body.mean
xy_product <- y_deviation * x_deviation

Очевидно, что величины (y_{i} — bar{y}) положительны для каждой точки в 1 и 2 квадрантах и отрицательны в 3 и 4 квадрантах. Аналогично (x_{i} — bar{x}) положительны в 1 и 4 квадрантах и отрицательны в 2 и 3 квадрантах. Если линейная связь между y и x отрицательна (с увеличением x уменьшается y), то больше точек в 2 и 4 квадрантах. И наоборот, если связь положительна (с увеличением x увеличивается y), то больше точек в 1 и 3 квадрантах. Таким образом, знак следующего выражения:

cov(y,x) = frac{sum_{i=1}^{n}(y_{i} — bar{y})(x_{i} — bar{x})}{n-1}

которое известно как коэффициент ковариации между y и x, указывает направление линейного отношения между y и x. Если cov(y,x) gt 0, то связь положительная, иначе отрицательная.

> sum(xy_product) / (nrow(animals.log) - 1)
[1] 7.051974

# Или с помощью функции cov()
> cov(x=animals.log$body, y=animals.log$brain)
[1] 7.051974

К сожалению, cov(Y,X) не может нам много сказать о силе связи, потому что она зависима от изменений в шкале измерения объектов. Для того, чтобы избежать этого недостатка ковариации необходимо нормализовать (стандартизировать) данные прежде чем вычислять ковариацию. Для нормализации Y сначала вычитается выборочное среднее из каждого наблюдения и делится на стандартное отклонение:

z_{i} = frac{y_{i} — bar{y}}{s_{y}}

где

s_{y} = sqrt{frac{sum_{i=1}^{n}(y_{i} — bar{y})^2}{n-1}}

называется выборочным среднеквадратическим отклонением Y. Легко показать, что нормализованная переменная Z имеет математическое ожидание 0 и стандартное отклонение 1. Аналогично нормализуется X. Ковариация между нормализованными X и Y называется коэффициентом корреляции между X и Y и вычисляется по следующей формуле:

cor(Y, X) = frac{1}{n-1}sum_{i=1}^{n}(frac{y_{i} — bar{y}}{s_{y}})(frac{x_{i} — bar{x}}{s_{x}}) = frac{cov(Y,X)}{s_{y}s_{x}} = frac{sum{(y_{i} — bar{y})(x_{i} — bar{x})}}{sqrt{sum{(y_{i} — bar{y})^2(x_{i} — bar{x})^2}}}

> cov(x=animals.log$body, y=animals.log$brain) / (sd(animals.log$body) * sd(animals.log$brain))
[1] 0.7794935

# Или с помощью функции cor()
> cor(animals.log$body, animals.log$brain)
[1] 0.7794935
par(mfrow=c(1, 2))
animals.log.cor <- cor(animals.log$body, animals.log$brain)
animals.cor <- cor(Animals$body, Animals$brain)
plot_animals(paste('r = ', round(animals.cor, 3)))
plot_animals(paste('r = ', round(animals.log.cor, 3)))
par(mfrow=c(1, 1))

Таким образом, Cor(Y, X) может быть интерпретирована как ковариация между стандартизированными переменными или отношение ковариации к стандартным отклонениям двух переменных. Коэффициент корреляции является симметричным, то есть: Cor(Y,X) = Cor(X,Y).

В отличии от Cov(Y,X), Cor(Y,X) инвариантна (не зависит от шкалы измеряемых величин) и удовлетворяет условию -1 leq Cor(Y,X) leq 1.

Итак, убедившись, что между y и x сущетвует линейная взаимосвязь необходимо построить модель, которая наилучшим образом позволила бы эту взаимосвязь описать, но какую модель следует выбрать?

plot_animals()
abline(lm(animals.log$brain ~ animals.log$body), col="darkred")
abline(a = 2.0, b = 0.8, col='darkred')
abline(a = 3.0, b = 0.3, col='darkred')

Мы рассмотрим разные линейные модели и начнем с простой линейной регрессии.

Todo

Рассказать про коэффициент корреляции, ближе к -1, 1 и 0. 0 не значит отсутствие зависимости, просто она может не быть линейной, например, Y = 50 — x^2. Рассказать про Anscombe датасет:

data(anscombe)
summary(anscombe)
cor(anscombe)

Простая линейная регрессия

Простая линейная регрессия (least squares regression line) это метод восстановления зависимости между двумя переменными: предсказываемой переменной Y и предиктором X:

Y = beta_{0} + beta_{1} X + epsilon

Где beta_{0} это пересечение (показывает значение Y при X = 0), beta_{1} угол наколна прямой, а epsilon случайная ошибка. Параметры beta_{0} и beta_{1} называются регрессионными коэффициентами.

Метод наименьших квадратов позволяет вычислить параметры beta_{0} и beta_{1} минимизируя сумму квадратов ошибок (вертикальных дистанций, см. рисунок ниже). Ошибки могут быть найдены по следующей формуле:

epsilon_{i} = y_{i} — beta_{0} — beta_{1}times x_{i}, i=1,2,…,n

Сумма квадратов ошибок (sum of squares errors) может быть записана как:

SSE(beta_{0}, beta_{1}) = sum_{i=1}^{n} epsilon_{i}^2 = sum_{i=1}^{n}(y_{i} — beta_{0} — beta_{1} x_{i})^2 rightarrow min

Note

В разной литературе можно встретить разные обозначения: SSE (sum of squared errors of prediction) == RSS (residual sum of squares) == SSR (sum of squared residuals).

Наша задача найти такие значения hat{beta_{0}} и hat{beta_{1}}, которые минимизируют SSE(beta_{0}, beta_{1}):

library(manipulate)

plotMininizeError <- function(beta0, beta1) {
    error <- sum((beta0 + beta1 * animals.log$body - animals.log$brain)^2)
    plot(animals.log$body, animals.log$brain,
         xlab='log(body weight)', ylab='log(brain weight)',
         main=paste(round(error, 3)))
    abline(a = beta0, b = beta1, col='darkred')
}

manipulate(plotMininizeError(b0, b1),
           b0=slider(0, 3, step=0.1),
           b1=slider(0, 2, step = 0.1))

Чтобы найти beta_{0}, beta_{1}, которые минимизируют SSE, необходимо решить:

triangledown SSE = left (frac{partial SSE}{partial beta_0}, frac{partial SSE}{partial beta_1}right ) = (0,0)

begin{aligned}
frac{partial SSE}{partial beta_0} =& -2sum_{i=1}^{n}(y_{i} — (beta_{0} + beta_{1} x_{i})) = 0 & Rightarrow &
&nbeta_0& + Big(sum_{i=1}^{n}x_iBig)beta_1 = sum_{i=1}^{n}y_i \
frac{partial SSE}{partial beta_1} =& -2sum_{i=1}^{n}(y_{i} — (beta_{0} + beta_{1} x_{i}))x_{i} = 0 & Rightarrow &
&Big(sum_{i=1}^{n}x_iBig)beta_0& + Big(sum_{i=1}^{n}x_i^2Big)beta_1 = sum_{i=1}^{n}x_iy_i
end{aligned}

Откуда получим:

beta_1 = frac{nsum x_iy_i — (sum x_i)(sum y_i)}{n(sum x_i^2) — (sum x_i)^2} = frac{sum (x_i-bar{x})(y_i-bar{y})}{sum (x_i-bar{x})} = rfrac{s_y}{s_x} qquad
beta_0 = bar{y} — beta_1bar{x}

Note

Линейная регрессия (или как ее еще иногда называют метод наименьших квадратов) была впервые описана Гауссом в 1795 году (Гауссу было всего 18 лет) и позже опубликована в его знаменитой книге «Theoria Motus Corporum Coelestium in sectionibus conicis solem ambientium» (Теория движения небесных тел), где он изложил методы вычисления планетных орбит.

В случае простой линейной регрессии коэффициент beta_{1} можно вычислить и более простым способом. Угол наклона прямой считается как frac{dy}{dx}, откуда получим:

begin{aligned}
beta_1 = frac{dy}{dx} rightarrow \
beta_1 = frac{sum_{i=1}^{n}(y_{i} — bar{y})}{sum_{i=1}^{n}(x_{i} — bar{x})} | times sum_{i=1}^{n}(x — bar{x}) rightarrow \
beta_1 = frac{sum_{i=1}^{n}(y_{i} — bar{y})(x_{i} — bar{x})}{sum_{i=1}^{n}(x_{i} — bar{x})^2} rightarrow … \
beta_1 = rfrac{s_{y}}{s_{x}}
end{aligned}

> (beta1 <- cor(animals.log$body, animals.log$brain) * sd(animals.log$brain) / sd(animals.log$body))
[1] 0.4959947
> (beta0 <- mean(animals.log$brain) - beta1 * mean(animals.log$body))
[1] 2.554898
> plot_animals(animals.log, xlab='log(body weight)', ylab='log(brain weight)',
               main = 'Brain and Body Weights for 28 Species')
> abline(a = beta0, b = beta1, col="darkred")

В R такую модель можно построить с помощью функции lm:

> model <- lm(formula = brain ~ body, data = animals.log)
> summary(model)

Call:
lm(formula = brain ~ body, data = animals.log)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2890 -0.6763  0.3316  0.8646  2.5835 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.55490    0.41314   6.184 1.53e-06 ***
body         0.49599    0.07817   6.345 1.02e-06 ***
---
Signif. codes:  0***0.001**0.01*0.05.’ 0.1 ‘ ’ 1

Residual standard error: 1.532 on 26 degrees of freedom
Multiple R-squared:  0.6076,    Adjusted R-squared:  0.5925 
F-statistic: 40.26 on 1 and 26 DF,  p-value: 1.017e-06

y ~ x называется формулой. Символ ~ в формуле следует читать как «описывается переменной». К объяснению вывода результатов команды summary мы еще вернемся. А пока ответим на вопрос «Как с помощью полученной модели делать прогнозы?».

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

hat{y} = hat{beta_0} + hat{beta_1}x

Подставив значения для оцененных коэффициентов hat{beta_0} и hat{beta_1} получим:

Теперь мы можем использовать эту модель для прогнозов.

library(manipulate)
plotPrediction <-function(newx) {
    model <- lm(brain~body, data=animals.log)
    plot(animals.log$body, animals.log$brain, 
         xlab='log(body weight)', ylab='log(brain weight)',
         main=paste('y = ', round(predict(model, newdata=data.frame(body=newx)), 3)))
    abline(model, lwd=2, col='dark blue')
    points(c(newx), c(min(animals.log$brain)-0.35), lwd=3, col='darkred')
    points(c(newx, newx), c(min(animals.log$brain)-0.35,
           model$coefficients[1] + newx*model$coefficients[2]),
           type='l', lwd=2, col="darkred", lty=2)
    points(c(min(animals.log$body)-0.4, newx),
           c(model$coefficients[1] + newx*model$coefficients[2],
           model$coefficients[1] + newx*model$coefficients[2]),
           type='l', lwd=2, col="darkred", lty=2)
    points(c(min(animals.log$body)-0.6),
           c(model$coefficients[1] + newx*model$coefficients[2]),
           lwd=3, col='darkred')
}

manipulate(plotPrediction(new_x), new_x = slider(min(animals.log$body), max(animals.log$body)))

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

begin{aligned}
ln(y) & = 2.555 + 0.496 times ln(x) + epsilon \
e^{ln(y)} & = e^{2.555 + 0.496 times ln(x) + epsilon} \
y & = e^{2.555}e^{ln(x)^{0.496}}e^{epsilon} \
y & = e^{2.555 + epsilon} times x^{0.496}
end{aligned}

Оценка полученной модели

Коэффициент детерминации (R-square) — это доля дисперсии зависимой переменной, объясняемая рассматриваемой моделью. Более точно — это единица минус доля необъяснённой дисперсии (дисперсии случайной ошибки модели, или условной по признакам дисперсии зависимой переменной) в дисперсии зависимой переменной:

R^{2} = 1 — frac{sum_{i=1}^{n}(y_i-hat{y_{i}})^2}{sum_{i=1}^{n}(y_i-bar{y})^2}

model <- lm(animals.log$brain ~ animals.log$body)
(r2 <- 1 - sum((animals.log$brain - predict(model))^2) / sum((animals.log$brain - mean(animals.log$brain))^2))
[1] 0.6076101

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

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

R{_{adj}}^{2} = 1 — frac{(1-R^2)(N-1)}{N-p-1}

(r2_adj <- 1 - ((1-r2)*(nrow(animals.log) - 1)) / (nrow(animals.log) - 1 - 1))
[1] 0.5925182

Также для оценки модели используют среднеквадратичную ошибку, которая характеризует отклонение реальных данных от линии регрессии и измеряется в тех же единицах, что и зависимая переменная (y):

RMSE = sqrt{frac{sum_{i=1}^{m}(hat{y_i} — y_i)^2}{m}}

Квартет Энскомба

Связь МНК и ММП

Рассмотрим принцип максимального правдоподобия (maximum likelihood estimation) на примере задачи подбрасывания монетки, а затем применим его к задаче линейной регрессии. Пусть у нас имеется:

  • y — число решек при подбрасывании монетки n раз;
  • n — число бросков;
  • p — вероятность выпадания решки (неизвестна);
  • hat{p} — оценка вероятности выпадания решки.

Предположим, что мы не знаем истинную вероятность выпадания решки (0.5), поэтому попробуем ее оценить. По интуиции:

hat{p} = frac{y}{n} = frac{# textit{решек}}{# textit{бросков}}

Например, если выпало 3 решки при 5 бросках, то:

hat{p} = frac{y}{n} = frac{3}{5} = 0.6

Теперь давайте попробуем подтвердить нашу интуицию математически. Эксперимент с подбрасыванием монетки может быть описан с помощью Биномиального распределения:

P(Y = y|n, p) = begin{pmatrix} n \ y end{pmatrix}p^y(1 — p)^{n — y}

Читаем как: «Вероятность того, что решка выпадет y раз при заданной вероятности p и числе бросков n».

Для y=3 и n=5 получим:

P(Y = 3|5, p) = begin{pmatrix} 5 \ 3 end{pmatrix}p^3(1 — p)^{5 — 3}

Но мы все еще не знаем p. Принцип максимального правдоподобия заключается в том, чтобы найти такое значение p для имеющихся y и n, при котором мы получим максимальное значение для P.

Итак, зафиксируем y и n, тогда функция правдоподобия (likelihood) будет выглядеть следующим образом:

L(p|y, n) = begin{pmatrix} n \ y end{pmatrix}p^y(1 — p)^{n — y} rightarrow max

Для наглядности будем перебирать значения p in [0;1] с шагом 0.1 и для каждого значения вычислим функцию правдоподобия:

library(graphics)
probs <- seq(0,1,0.1)
likelihoods <- dbinom(3, size=5, prob=probs)
plot(probs, likelihoods, type='l', xlab = 'probability', ylab = 'likelihood')
points(probs, likelihoods, col = 'royal blue')
segments(probs, rep(0, length(probs)), probs, likelihoods, lty=3, col='royal blue')

Из графика хорошо видно, что мы достигаем максимума при p = 0.6, но мы не будем каждый раз перебирать параметр p.

Иногда для упрощения расчетов берут логарифм, таким образом, максимизируя log-likelihood:

log ell (p|y, n) = log begin{pmatrix} n \ y end{pmatrix} + ylog p + (n — y)log (1 — p)

И наконец продифференцируем функцию относительно параметра p и приравняем ее нулю:

frac{partial log ell}{partial p} = frac{y}{p} — frac{n — y}{1 — p} = 0 Rightarrow y — yp — np + yp = 0 Rightarrow hat{p} = frac{y}{n}

Теперь кратко рассмотрим метод максимального правдоподобия применительно к линейной регрессии:

Y_i = beta_0 + beta_1x_i + epsilon

Предположим, что остатки имеют нормальное распределение epsilon_i sim mathcal{N}(0, sigma^2), тогда:

Y_i|X_i sim mathcal{N}(beta_0 + beta_1 x_i, sigma^2)

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

f_X(x; mu, sigma^2) = frac{1}{sqrt{2pisigma^2}}exp -frac{(x — mu)^2}{2sigma^2}

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

f_{Y_i}(y_i; x_i, beta_0, beta_1, sigma^2) = frac{1}{sqrt{2pisigma^2}}exp -frac{(y_i — beta_0 — beta_1 x_i)^2}{2sigma^2}

Для n точек (x_i, y_i) мы можем записать функцию правдоподобия от трех параметров beta_0, beta_1, sigma^2 следующим образом:

L(beta_0, beta_1, sigma^2) = prod_{i = 1}^{n}frac{1}{sqrt{2pisigma^2}}exp -frac{(y_i — beta_0 — beta_1 x_i)^2}{2sigma^2}

Наша задача найти такие параметры beta_0, beta_1 и sigma^2, которые максимизируют L(beta_0, beta_1, sigma^2). Для простоты вычислений будем считать log-likelihood:

Note

По интуиции должно быть понятно, что log ell rightarrow max это тоже самое, что и SSE(beta_0, beta_1) rightarrow min.

begin{align}
log ell =& -sum_{i = 1}^{n}Big(frac{1}{2}log2pi + frac{1}{2}logsigma^2 + frac{1}{2sigma^2}(y_i — beta_0 — beta_1 x_i)^2Big) \
=& -frac{n}{2}log2pi -frac{n}{2}logsigma^2 — frac{1}{2sigma^2}sum_{i=1}^{n}(y_i — beta_0 — beta_1 x_i)^2
end{align}

Как и в задаче с подбрасыванием монетки возьмем частные проиводные по параметрам beta_0, beta_1 и sigma^2 и затем приравняем нулю:

frac{partial log ell}{partial beta_0} = frac{1}{sigma^2}sum_{i=1}^{n}(y_i — beta_0 — beta_1 x_i) Rightarrow sum_{i=1}^{n}(y_i — beta_0 — beta_1 x_i) = 0

frac{partial log ell}{partial beta_1} = frac{1}{sigma^2}sum_{i=1}^{n}(y_i — beta_0 — beta_1 x_i)x_i Rightarrow sum_{i=1}^{n}(y_i — beta_0 — beta_1 x_i)x_i = 0

frac{partial log ell}{partial sigma^2} = -frac{n}{2 sigma^2} + frac{1}{2 (sigma^2)^2}sum_{i=1}^{n}(y_i — beta_0 — beta_1 x_i)^2 Rightarrow -frac{n}{2 sigma^2} + frac{1}{2 (sigma^2)^2}sum_{i=1}^{n}(y_i — beta_0 — beta_1 x_i)^2 = 0

Откуда получим:

hat{beta_0} = bar{y} — beta_1bar{x}

hat{beta_1} = rfrac{s_y}{s_x}

hat{sigma}^2 = frac{1}{n}sum_{i=1}^{n}(y_i — hat{y_i})^2

Многомерная линейная регрессия

y_i = beta_1 + beta_2x_{2i} + … + beta_kx_{ki} + epsilon_i,quad (i=1,…,n)

Запишем в матричной форме:

y = begin{pmatrix} y_1 \ … \ y_n end{pmatrix},quad
X = begin{pmatrix} 1 & x_{21} & … & x_{k1} \ … & … & … & … \ 1 & x_{2n} & … & x_{kn} end{pmatrix},quad
beta = begin{pmatrix} beta_1 \ … \ beta_k end{pmatrix},quad
epsilon = begin{pmatrix} epsilon_1 \ … \ epsilon_n end{pmatrix}

Опять хотим минимизировать квадрат ошибки:

begin{align*}
L = sum_{i=1}^{n}epsilon_{i}^2 = epsilon’epsilon &= (y — Xbeta)'(y-Xbeta) \
&= y’y — y’Xbeta — X’beta’y + beta’X’b rightarrow min
end{align*}

Уже по известной схеме:

frac{partial L}{partial beta} = -2X’y + 2X’Xbeta = 0

Откуда получим:

X’Xb = X’y Rightarrow beta = (X’X)^{-1}X’y

Последнее уравнение называется нормальным уравнением (normal equation).

n <- length(y)
X <- cbind(x, rep(1,n))
B <- solve(t(X)%*%X)%*%t(X)%*%y
x.test <- seq(1,5, 0.25)
y.hat <- B[2]+B[1]*x.test
lines(x.test, y.hat,col="red")

Робастность регрессии

par(mfrow=c(1,2))
plot(brain ~ body, data=animals.log,
     xlab='log(body weight)', ylab='log(brain weight)')
points(animals.log[animals.log$body > 8 & animals.log$brain < 6,],
       pch=16, col="darkred")
model <- lm(brain ~ body, data = animals.log)
abline(model, lwd=2, col='dark blue')

plot(brain ~ body, data=animals.log,
     xlab='log(body weight)', ylab='log(brain weight)')
points(animals.log[animals.log$body > 8 & animals.log$brain < 6,],
       pch=16, col="darkred")
model <- rlm(brain ~ body, data=animals.log)
abline(model, lwd=2, col='dark blue')
par(mfrow=c(1,1))

Полиномиальная регрессия, ridge-регрессия и LASSO-регрессия

par(mfrow=c(1,2))
set.seed(123456)
x <- seq(-2, 5, length.out=20)
y <- sin(x) + runif(length(x), -0.25, 0.25)
fit <- lm(y ~ x)
plot(y ~ x, lwd=2, main="Ordinary Least Squares")
abline(fit, lwd=2, col='dark blue')
plot(resid(fit) ~ fitted(fit), lwd=2, ylab="Residuals", xlab="Fitted", main="Residuals vs Fitted values")
abline(0, 0, lwd=2, col='dark blue')
par(mfrow=c(1,1))

par(mfrow=c(3,3))

poly.1  <- lm(y ~ poly(x, 1))
poly.7  <- lm(y ~ poly(x, 7))
poly.14 <- lm(y ~ poly(x, 14))

plot(y ~ x, lwd = 2, main="")
curve(predict(poly.1,  data.frame(x=x)), add=TRUE, col="royalblue")
plot(y ~ x, lwd = 2, main="")
curve(predict(poly.7,  data.frame(x=x)), add=TRUE, col="darkgreen")
plot(y ~ x, lwd = 2, main="")
curve(predict(poly.14, data.frame(x=x)), add=TRUE, col="darkred")

plot(y ~ x, lwd = 2, main="")
curve(predict(poly.1,  data.frame(x=x)), add=TRUE, col="royalblue")
points(x1, y1, lwd=2, pch=2)
plot(y ~ x, lwd = 2, main="")
curve(predict(poly.7,  data.frame(x=x)), add=TRUE, col="darkgreen")
points(x1, y1, lwd=2, pch=2)
plot(y ~ x, lwd = 2, main="")
curve(predict(poly.14, data.frame(x=x)), add=TRUE, col="darkred")
points(x1, y1, lwd=2, pch=2)

plot(resid(poly.1) ~ fitted(poly.1), lwd=2, ylab="Residuals", xlab="Fitted")
abline(0, 0, lwd=2, col='dark blue')
plot(resid(poly.7) ~ fitted(poly.7), lwd=2, ylab="Residuals", xlab="Fitted")
abline(0, 0, lwd=2, col='dark blue')
plot(resid(poly.14) ~ fitted(poly.14), lwd=2, ylab="Residuals", xlab="Fitted")
abline(0, 0, lwd=2, col='dark blue')

par(mfrow=c(1,1))

Простая модель очень плохо описывает наши данные и это называется underfitting, полиномиальная модель же наоборот, прошла через каждую точку и в точности описала все данные, это называется overfitting.

poly.2 <- lm(y ~ poly(x, 2))
poly.3 <- lm(y ~ poly(x, 3))
poly.4 <- lm(y ~ poly(x, 4))
poly.5 <- lm(y ~ poly(x, 5))
poly.6 <- lm(y ~ poly(x, 6))
poly.8 <- lm(y ~ poly(x, 8))
poly.9 <- lm(y ~ poly(x, 9))
poly.10 <- lm(y ~ poly(x, 10))
poly.11 <- lm(y ~ poly(x, 11))
poly.12 <- lm(y ~ poly(x, 12))
poly.13 <- lm(y ~ poly(x, 13))
anova(poly.1, poly.2, poly.3, poly.4, poly.5, poly.6, poly.7, poly.8, poly.9, 
      poly.10, poly.11, poly.12)
Analysis of Variance Table

Model  1: y ~ poly(x, 1)
Model  2: y ~ poly(x, 2)
Model  3: y ~ poly(x, 3)
Model  4: y ~ poly(x, 4)
Model  5: y ~ poly(x, 5)
Model  6: y ~ poly(x, 6)
Model  7: y ~ poly(x, 7)
Model  8: y ~ poly(x, 8)
Model  9: y ~ poly(x, 9)
Model 10: y ~ poly(x, 10)
Model 11: y ~ poly(x, 11)
Model 12: y ~ poly(x, 12)
   Res.Df    RSS Df Sum of Sq       F  Pr(>F)  
   1      13 6.7847                               
   2      12 1.7598  1    5.0249 68.7312 0.01424 *
   3      11 1.6804  1    0.0794  1.0859 0.40679  
   4      10 0.2022  1    1.4782 20.2190 0.04607 *
   5       9 0.1992  1    0.0030  0.0416 0.85728  
   6       8 0.1548  1    0.0444  0.6075 0.51732  
   7       7 0.1542  1    0.0006  0.0079 0.93723  
   8       6 0.1542  1    0.0000  0.0004 0.98568  
   9       5 0.1473  1    0.0069  0.0937 0.78845  
   10      4 0.1470  1    0.0003  0.0042 0.95421  
   11      3 0.1469  1    0.0001  0.0017 0.97091  
   12      2 0.1462  1    0.0007  0.0090 0.93296  
   ---
   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Ridge-регрессия

Note

Регуляризация Тихонова.

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

Ранее мы минимизировали:

SSE(beta) = sum_{i=1}^{n}(y_{i} — sum_{j=1}^{p}beta_{j} x_{j})^2

А теперь будем минимизировать:

SSE(beta) = sum_{i=1}^{n}(y_{i} — sum_{j=1}^{p}beta_{j} x_{j})^2 + lambdasum_{j=1}^{p}beta_{j}^2

где lambda — параметр регуляризации, который надо как-то выбрать.

L = (y — Xbeta)^{T}(y-Xbeta) + lambdabeta^{T}beta

Если взять производную, то получим:

beta = (X^{T}X + lambda I)^{-1}X^{T}y

x <- seq(-2, 5, length.out=15)
y <- sin(x) + runif(length(x), -0.25, 0.25)
plot(y ~ x, lwd = 2)

x1 <- seq(-2, 5, length.out = 20)
y <- sin(x1) + runif(length(x1), -0.25, 0.25)
points(x1, y1, col="darkred", lwd=2, pch=17)

cv.ridge <- cv.glmnet(poly(x, 7), y, alpha=0, family='gaussian')
points(x1, predict(cv.ridge, poly(x1, 7), s=cv.ridge$lambda.min), type="l", col="royalblue", lwd=2)

LASSO-регрессия

SSE(beta) = sum_{i=1}^{n}(y_{i} — sum_{j=1}^{p}beta_{j} x_{j})^2 + lambdasum_{j=1}^{p}beta_{j}^2

library(lars)
lasso.model <- lars(x = poly(x, 14), y = y, type = "lasso")
plot(lasso.model)

x1 <- seq(-2, 5, length.out = 20)
y1 <- sin(x1) + runif(20, -0.25, 0.25)
y.pred.lasso <- predict(lasso.model,
                        poly(x1, 14),
                        s=best.fraction,
                        type="fit",
                        mode="fraction")
plot(y ~ x, lwd = 2, main="")
points(x1, y.pred.lasso$fit, type="l", col="royalblue", lwd=2)
points(x1, y1, col="darkred", lwd=2, pch=17)


Последнее обновление: 4 мая 2021 г.

0 0 голоса
Рейтинг статьи
Подписаться
Уведомить о
guest

0 комментариев
Старые
Новые Популярные
Межтекстовые Отзывы
Посмотреть все комментарии

А вот еще интересные материалы:

  • Яшка сломя голову остановился исправьте ошибки
  • Ясность цели позволяет целеустремленно добиваться намеченного исправьте ошибки
  • Ясность цели позволяет целеустремленно добиваться намеченного где ошибка
  • Лексус код ошибки p0352
  • Лефт фо дед 2 ошибка стим