Как найти функцию методом наименьших квадратов

Метод наименьших квадратов: формулы, код и применение

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

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

Обычно объяснение метода представляется в виде:

1. перегруженного математическими абстракциями «классического описания»;

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

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

Триада

Данные проблемы возникают из-за того, что при описании МНК не описывается путь перехода от математического описания метода к его программной реализации. Данная проблема была давно преодолена одним из отцов-основателей математического моделирования в нашей стране академиком Самарским А.А., который предложил универсальную методологию вычислительного эксперимента в виде триады «модель-алгоритм-программа».

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

1. описать математическую модель объекта/процесса, анализируемые параметры которой формализуются в виде некой функциональной зависимости;

2. синтезировать алгоритм оценки искомых параметров модели на основе МНК;

3. выполнить программную реализацию полученного алгоритма.

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

В качестве анализируемых данных будем рассматривать n значений температуры, измеренной в градусах, получаемые с использованием термометра и зависящие от времени. Эти данные представим в виде таблицы, где строки – это объекты, первый столбец – набор параметров (j=1), второй столбец – набор измеренных ответов:

т.е. случайная величина, распределённая по закону Гаусса с фиксированным значением среднеквадратического отклонения (СКО), нулевым математическим ожиданием (МО) и с автокорреляционной функцией в виде дельта-функции.

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

import random
import matplotlib.pyplot as plt
import numpy as np

def datasets_make_regression(coef, data_size, noise_sigma, random_state):
    x = np.arange(0, data_size, 1.)
    mu = 0.0
    random.seed(random_state)
    noise = np.empty((data_size, 1))
    y = np.empty((data_size, 1))

    for i in range(data_size):
        noise[i] = random.gauss(mu, noise_sigma)
        y[i] = coef[0] + coef[1]*x[i] + noise[i]

    return x, y

coef_true = [34.2, 2.] # весовые коэффициенты
data_size = 200 # размер генерируемого набора данных
noise_sigma = 10 # СКО шума в данных
random_state = 42
x_scale, y_estimate = datasets_make_regression(coef_true, data_size, noise_sigma, random_state)

plt.plot(x_scale, y_estimate, 'o') 
plt.xlabel('x (порядковый номер измерения)', fontsize=14)
plt.ylabel('y (оценка температуры)', fontsize=14)

Сгенерированный набор данных будет иметь вид

Математическая модель

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

Первое описание математической модели. Самое общее.

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

Второе описание математической модели. Общее описание линейной модели.

Поскольку по условию задачи мы рассматриваем линейную зависимость, на которую влияет шум, имеющий нормальное распределение с нулевым средним значением, то используем линейную относительно параметров модель, а метрикой выберем расчёт СКО ошибки

Третье описание математической модели. Модель конкретизирована под условия примера.

Поскольку по условию задачи мы рассматриваем линейную зависимость одного параметра (поскольку j=1, то далее эту зависимость при описании опустим), на оценку которого влияет шум, имеющий нормальное распределение, то используем линейную модель с двумя коэффициентами

Обсуждение моделей

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

Хочется отметить, что параметры математической модели процесса изменения температуры, описанной в выражении (2), в принципе, могут находиться не только методом наименьших квадратов, поскольку метрика с квадратичной функцией ошибки в явном виде не прописана. Такая общая форма записи приведена, поскольку не представляется возможным в кратком виде привести все возможные модификации МНК, у которого минимизация для различных условий может выполняться также различно.

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

Далее с использованием выражения (4) и МНК, синтезируем алгоритмы, оценивающие параметры модели.

Алгоритмы и реализующие их программы

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

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

применим правило дифференцирования сложной функции

применим правило дифференцирования функции представленной в виде суммы

применим правило дифференцирования функции представленной в виде суммы

Приравняем полученные производные к нулю и решим полученную систему уравнений

Раскроем скобки

Вынесем постоянные множители за скобки

Вынесем слагаемые с множителем «y» в правую часть уравнений

Поставим слагаемые с множителем «x» в левой части в порядке убывания степеней

Для решения полученной системы алгебраических уравнения представим её в матричной форме

Выразим вектор w* с искомыми весами выполнив умножение обеих частей равенства на обратную матрицу y

Полученное выражение (5) является решением системы уравнений и его можно уже использовать в качестве первого алгоритма оценки параметров модели. Ниже приведен код реализующий этот алгоритм.

def coefficient_reg_inv(x, y):
    size = len(x)
    # формируем и заполняем матрицу размерностью 2x2
    A = np.empty((2, 2))
    A[[0], [0]] = sum((x[i])**2 for i in range(0,size))
    A[[0], [1]] = sum(x)
    A[[1], [0]] = sum(x)
    A[[1], [1]] = size
    # находим обратную матрицу 
    A = np.linalg.inv(A)
    # формируем и заполняем матрицу размерностью 2x1
    C = np.empty((2, 1))
    C[0] = sum((x[i]*y[i]) for i in range(0,size))
    C[1] = sum((y[i]) for i in range(0,size))
    
    # умножаем матрицу на вектор
    ww = np.dot(A, C)
    return ww[1], ww[0]

[w0_1, w1_1] = coefficient_reg_analit(x_scale, y_estimate)
print(w0_1, w1_1)

Результат работы алгоритма:
[33.93193341] [2.01436546]

Выражение (5) можно упростить, выполнив аналитический расчет обратной матрицы.

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

Перепишем выражение (5) в виде

Упростим выражение раскрыв скобки, чтобы получить более компактную форму записи

Теперь решение системы уравнений имеет вид

Полученное выражение (6) можно использовать в качестве второго алгоритма оценки параметров модели. Ниже приведен код, реализующий этот алгоритм.

def coefficient_reg_inv_analit(x, y):
    size = len(x)
    # выполним расчет числителя первого элемента вектора
    numerator_w1 = size*sum(x[i]*y[i] for i in range(0,size)) - sum(x)*sum(y)
    # выполним расчет знаменателя (одинаковый для обоих элементов вектора)
    denominator = size*sum((x[i])**2 for i in range(0,size)) - (sum(x))**2
    # выполним расчет числителя второго элемента вектора
    numerator_w0 = -sum(x)*sum(x[i]*y[i] for i in range(0,size)) + sum((x[i])**2 for i in range(0,size))*sum(y)
    
    # расчет искомых коэффициентов
    w1 = numerator_w1/denominator
    w0 = numerator_w0/denominator
    return w0, w1

[w0_2, w1_2] = coefficient_reg_inv_analit(x_scale, y_estimate)
print(w0_2, w1_2)

Результат работы алгоритма:
[33.93193341] [2.01436546]

Видно, что полученный результат совпадает с результатом, полученным с использованием выражения (5).

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

раскроем скобки:

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

С учётом введённых переменных искомый вектор w* примет вид

Таким образом итоговое выражение представим в виде:

Полученное выражение (7) можно использовать в качестве третьего алгоритмаоценки параметров модели. Ниже приведен код, реализующий этот алгоритм.

def coefficient_reg_stat(x, y):
    size = len(x)
    avg_x = sum(x)/len(x) # оценка МО величины x
    avg_y = sum(y)/len(y) # оценка МО величины y
    # оценка МО величины x*y
    avg_xy = sum(x[i]*y[i] for i in range(0,size))/size
    # оценка СКО величины x 
    std_x = (sum((x[i] - avg_x)**2 for i in range(0,size))/size)**0.5
    # оценка СКО величины y
    std_y = (sum((y[i] - avg_y)**2 for i in range(0,size))/size)**0.5
    # оценка коэффициента корреляции величин x и y
    corr_xy = (avg_xy - avg_x*avg_y)/(std_x*std_y)
    
    # расчет искомых коэффициентов
    w1 = corr_xy*std_y/std_x
    w0 = avg_y - avg_x*w1
    return w0, w1

[w0_3, w1_3] = coefficient_reg_stat(x_scale, y_estimate)
print(w0_3, w1_3)

Результат работы алгоритма:
[33.93193341] [2.01436546]

Сравним полученные значения с результатом работы процедуры библиотеки sklearn.

from sklearn.linear_model import LinearRegression
# преобразование размерности массива x_scale для корректной работы model.fit
x_scale = x_scale.reshape((-1,1)) 
model = LinearRegression()
model.fit(x_scale, y_estimate)
print(model.intercept_, model.coef_)

Результат работы алгоритма:
[33.93193341] [2.01436546]

Видно, что полученный результат полностью совпадает с результатами, полученными ранее.

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

def predict(w0, w1, x_scale):
    y_pred = [w0 + val*w1 for val in x_scale]
    return y_pred

y_predict = predict(w0_1, w1_1, x_scale)

plt.plot(x_scale, y_estimate, 'o', label = 'Истинные значения')
plt.plot(x_scale, y_predict, '*', label = 'Расчетные значения')
plt.legend(loc = 'best', fontsize=12) 
plt.xlabel('x (порядковый номер измерения)', fontsize=14)
plt.ylabel('y (оценка температуры)', fontsize=14)

Обсуждение алгоритмов

Таким образом удалось выполнить оценку параметров модели, которая описывает процесс изменения температуры в печи с минимальной ошибкой. Выведенные выражения (5) – (7) и запрограммированные по ним три алгоритма, реализующие МНК при обработке одного набора измеренных температур, хотя и отличаются видом конечных выражений, дают одинаковые оценки. Последовательный вывод выражений, выполненный в статье, показывает, что эти выражения по сути «делают» одно и тоже, однако, позволяют давать различные интерпретации. При этом хочется отметить, что в алгоритмах (6) – (7) удалось уйти от процедуры обращения матрицы, которая при обработке реальных данных может выполняться нестабильно.

Ещё один пример. Анализируем остатки.

• количество выбросов в выборке равно 10;

• длительность импульса равна 1-му измерению;

• амплитуда импульса равна 256;

• позиция первого импульса соответствует 10-му измерению;

• период следования импульсов – каждые 20 измерений.

ni = 10 # количество выбросов
ind_impuls = np.arange(ni, data_size, 20) # индексы выбросов
y_estimate_imp = y_estimate.copy() # выборка с выбросами
for i in range(0, ni):
    y_estimate_imp[ind_impuls[i]] += 256

[w0_imp, w1_imp] = coefficient_reg_stat(x_scale, y_estimate_imp)
y_pred_imp = predict(w0_imp, w1_imp, x_scale)
plt.plot(x_scale, y_estimate_imp, 'o', label = 'Истинные значения')
plt.plot(x_scale, y_pred_imp, '*', label = 'Расчетные значения')
plt.legend(loc = 'best', fontsize=12) 
plt.xlabel('x (порядковый номер измерения)', fontsize=14)
plt.ylabel('y (оценка температуры)', fontsize=14) 

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

Ниже приведен код, который выполняет этот расчет (8) и визуализирует полученный результат (рисунок 4).

SqErr = (y_pred_imp - y_estimate_imp)**2

plt.plot(x_scale, SqErr, 'o') 
plt.xlabel('x (порядковый номер измерения)', fontsize=14)
plt.ylabel('Квадрат ошибки', fontsize=14)

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

def censor_data(SqErr, nCensor):
    # индексы отсортированного во возрастанию массива с квадратами ошибок
    I = np.argsort(SqErr[:,0])
    ind_imp = I[-nCensor:]
    ind_imp = ind_imp[::-1] # разворот индексов массива
    w0 = np.empty((nCensor, 1))
    w1 = np.empty((nCensor, 1))

    for i in range(0,nCensor):
        # цензурирование данных
        x_scale_cens = np.delete(x_scale, ind_imp[0:i], 0)
        y_estimate_imp_cens = np.delete(y_estimate_imp, ind_imp[0:i], 0)
        # расчёт параметров модельной прямой
        w0[i], w1[i] = coefficient_reg_stat(x_scale_cens, y_estimate_imp_cens)
        y_pred2_cens = predict(w0[i], w1[i], x_scale_cens) 
    return w0, w1 

nCensor = 20 # количество отбрасываемых выбросов
[w0_с, w1_с] = censor_data(SqErr, nCensor)

plt.plot(coef_true[0], coef_true[1], 'o', label = 'Истинные значения')
plt.plot(w0_с, w1_с, '-*', label = 'Расчетные значения')
plt.legend(loc = 'best', fontsize=12)
plt.xlim((30, 50))
plt.ylim((1.85, 2.15))
plt.xlabel('w0', fontsize=18)
plt.ylabel('w1', fontsize=18)

По итогам работы алгоритма построим полученные оценки параметров модели (рисунок 5) по мере цензурирования данных.

В результате работы рассматриваемого алгоритма видно (рисунок 5), что по мере цензурирования помеховых данных параметры модели приближаются к истинным (точка с оценкой параметров возле i=1 получена после удаления первой точки, точка с оценкой параметров возле i=20 получена после удаления двадцатой точки).

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

y_pred_censor = predict(w0_с[nCensor-1], w1_с[nCensor-1], x_scale)

plt.plot(x_scale, y_estimate_imp, 'o', label = 'Истинные значения')
plt.plot(x_scale, y_pred_censor, '*', label = 'Расчетные значения')
plt.legend(loc = 'best', fontsize=12)
plt.xlabel('X (порядковый номер измерения)', fontsize=14)
plt.ylabel('Y (оценка температуры)', fontsize=14)

Обсуждение алгоритма

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

Безусловно, что приведенный в статье алгоритм оценки параметров модели с использованием процедуры цензурирования импульсных помеховых сигналов (в том виде, в котором он приведен в статье) не является оптимальным и универсальным. Он приведен здесь для иллюстрации возможностей модификации алгоритма, реализующего классический МНК. А также для демонстрации ещё одной важной метрики при использовании МНК – оценка квадратов разности между измеренными и модельными данными (8), которая является источником множества различных модификаций метода.

Заключение

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

Начнем статью сразу с примера. У нас есть некие экспериментальные данные о значениях двух переменных – x и y. Занесем их в таблицу.

После выравнивания получим функцию следующего вида: g(x)=x+13+1.

Мы можем аппроксимировать эти данные с помощью линейной зависимости y=ax+b, вычислив соответствующие параметры. Для этого нам нужно будет применить так называемый метод наименьших квадратов. Также потребуется сделать чертеж, чтобы проверить, какая линия будет лучше выравнивать экспериментальные данные.

В чем именно заключается МНК (метод наименьших квадратов)

Главное, что нам нужно сделать, – это найти такие коэффициенты линейной зависимости, при которых значение функции двух переменных F(a, b)=∑i=1n(yi-(axi+b))2 будет наименьшим. Иначе говоря, при определенных значениях a и b сумма квадратов отклонений представленных данных от получившейся прямой будет иметь минимальное значение. В этом и состоит смысл метода наименьших квадратов. Все, что нам надо сделать для решения примера – это найти экстремум функции двух переменных.

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

Для того чтобы вывести формулы для вычисления коэффициентов, нужно составить и решить систему уравнений с двумя переменными. Для этого мы вычисляем частные производные выражения F(a, b)=∑i=1n(yi-(axi+b))2 по a и b и приравниваем их к 0.

δF(a, b)δa=0δF(a, b)δb=0⇔-2∑i=1n(yi-(axi+b))xi=0-2∑i=1n(yi-(axi+b))=0⇔a∑i=1nxi2+b∑i=1nxi=∑i=1nxiyia∑i=1nxi+∑i=1nb=∑i=1nyi⇔a∑i=1nxi2+b∑i=1nxi=∑i=1nxiyia∑i=1nxi+nb=∑i=1nyi

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

n∑i=1nxiyi-∑i=1nxi∑i=1nyin∑i=1n-∑i=1nxi2b=∑i=1nyi-a∑i=1nxin

Мы вычислили значения переменных, при который функция
F(a, b)=∑i=1n(yi-(axi+b))2 примет минимальное значение. В третьем пункте мы докажем, почему оно является именно таким.

Это и есть применение метода наименьших квадратов на практике. Его формула, которая применяется для поиска параметра a, включает в себя ∑i=1nxi, ∑i=1nyi, ∑i=1nxiyi, ∑i=1nxi2, а также параметр
n – им обозначено количество экспериментальных данных. Советуем вам вычислять каждую сумму отдельно. Значение коэффициента b вычисляется сразу после a.

Обратимся вновь к исходному примеру.

Как изобразить МНК на графике функций

Метод наименьших квадратов наглядно показан на графической иллюстрации. С помощью красной линии отмечена прямая g(x)=x+13+1, синей – y=0,165x+2,184. Исходные данные обозначены розовыми точками.

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

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

Доказательство метода МНК

Чтобы функция приняла минимальное значение при вычисленных a и b, нужно, чтобы в данной точке матрица квадратичной формы дифференциала функции вида F(a, b)=∑i=1n(yi-(axi+b))2 была положительно определенной. Покажем, как это должно выглядеть.

Запрос «МНК» перенаправляется сюда; см. также другие значения.

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

История[править | править код]

До начала XIX в. учёные не имели определённых правил для решения системы уравнений, в которой число неизвестных меньше, чем число уравнений; до этого времени употреблялись частные приёмы, зависевшие от вида уравнений и от остроумия вычислителей, и потому разные вычислители, исходя из тех же данных наблюдений, приходили к различным выводам. Гауссу (1795) принадлежит первое применение метода, а Лежандр (1805) независимо открыл и опубликовал его под современным названием (фр. Méthode des moindres quarrés)^[1]. Лаплас связал метод с теорией вероятностей, а американский математик Эдрейн^ru_en (1808) рассмотрел его теоретико-вероятностные приложения^[2]. Метод распространён и усовершенствован дальнейшими изысканиями Энке, Бесселя, Ганзена и других.

Работы А. А. Маркова в начале XX века позволили включить метод наименьших квадратов в теорию оценивания математической статистики, в которой он является важной и естественной частью. Усилиями Ю. Неймана, Ф. Дэвида, А. Эйткена, С. Рао было получено множество немаловажных результатов в этой области^[3].

Суть метода наименьших квадратов[править | править код]

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

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

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

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

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

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

Пример — система линейных уравнений[править | править код]

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

,

где прямоугольная матрица размера (то есть число строк матрицы A больше количества искомых переменных).

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

.

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

,

где  — псевдообратная матрица для .

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

Строгое обоснование и установление границ содержательной применимости метода даны А. А. Марковым и А. Н. Колмогоровым.

МНК в регрессионном анализе (аппроксимация данных)[править | править код]

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

.

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

,

где  — так называемые случайные ошибки модели.

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

,

где  — англ. Residual Sum of Squares^[4] определяется как:

.

В общем случае решение этой задачи может осуществляться численными методами оптимизации (минимизации). В этом случае говорят о нелинейном МНК (NLS или NLLS — англ. Non-Linear Least Squares). Во многих случаях можно получить аналитическое решение. Для решения задачи минимизации необходимо найти стационарные точки функции , продифференцировав её по неизвестным параметрам , приравняв производные к нулю и решив полученную систему уравнений:

.

МНК в случае линейной регрессии[править | править код]

Пусть регрессионная зависимость является линейной:

.

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

.

Тогда вектор оценок объясняемой переменной и вектор остатков регрессии будут равны

.

соответственно сумма квадратов остатков регрессии будет равна

.

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

.

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

где все суммы берутся по всем допустимым значениям .

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

Решение этой системы уравнений и дает общую формулу МНК-оценок для линейной модели:

.

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

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

.

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

Простейшие частные случаи[править | править код]

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

.

Отсюда несложно найти оценки коэффициентов:

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

.

Следовательно, формула оценки единственного коэффициента имеет вид

.

Случай полиномиальной модели[править | править код]

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

Статистические свойства МНК-оценок[править | править код]

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

1. математическое ожидание случайных ошибок равно нулю и
2. факторы и случайные ошибки — независимые случайные величины.

Первое условие для моделей с константой можно считать выполненным всегда, так как константа берёт на себя ненулевое математическое ожидание ошибок (поэтому модели с константой в общем случае предпочтительнее).

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

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

• Постоянная (одинаковая) дисперсия случайных ошибок во всех наблюдениях (отсутствие гетероскедастичности): .

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

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

Линейная модель, удовлетворяющая таким условиям, называется классической. МНК-оценки для классической линейной регрессии являются несмещёнными, состоятельными и наиболее эффективными оценками в классе всех линейных несмещённых оценок (в англоязычной литературе иногда употребляют аббревиатуру BLUE (Best Linear Unbiased Estimator) — наилучшая линейная несмещённая оценка; в отечественной литературе чаще приводится теорема Гаусса — Маркова). Как нетрудно показать, ковариационная матрица вектора оценок коэффициентов будет равна:

.

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

.

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

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

Обобщённый МНК[править | править код]

Метод наименьших квадратов допускает широкое обобщение. Вместо минимизации суммы квадратов остатков можно минимизировать некоторую положительно определённую квадратичную форму от вектора остатков , где  — некоторая симметрическая положительно определённая весовая матрица. Обычный МНК является частным случаем данного подхода, когда весовая матрица пропорциональна единичной матрице. Как известно, для симметрических матриц (или операторов) существует разложение . Следовательно, указанный функционал можно представить следующим образом: , то есть этот функционал можно представить как сумму квадратов некоторых преобразованных «остатков».
Таким образом, можно выделить класс методов наименьших квадратов — LS-методы (Least Squares).

Доказано (теорема Айткена), что для обобщённой линейной регрессионной модели (в которой на ковариационную матрицу случайных ошибок не налагается никаких ограничений) наиболее эффективными (в классе линейных несмещённых оценок) являются оценки т. н. обобщённого МНК (ОМНК, GLS — Generalized Least Squares) — LS-метода с весовой матрицей, равной обратной ковариационной матрице случайных ошибок: .

Можно показать, что формула ОМНК-оценок параметров линейной модели имеет вид

.

Ковариационная матрица этих оценок соответственно будет равна

.

Фактически сущность ОМНК заключается в определённом (линейном) преобразовании (P) исходных данных и применении обычного МНК к преобразованным данным. Цель этого преобразования — для преобразованных данных случайные ошибки уже удовлетворяют классическим предположениям.

Взвешенный МНК[править | править код]

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

См. также[править | править код]

• Обобщенный метод наименьших квадратов
• Двухшаговый метод наименьших квадратов
• Рекурсивный МНК
• Алгоритм Гаусса — Ньютона

Примечания[править | править код]

Литература[править | править код]

• Линник Ю. В. Метод наименьших квадратов и основы математико-статистической теории обработки наблюдений. — 2-е изд. — М., 1962. (математическая теория)
• Айвазян С. А. Прикладная статистика. Основы эконометрики. Том 2. — М.: Юнити-Дана, 2001. — 432 с. — ISBN 5-238-00305-6.
• Доугерти К. Введение в эконометрику: Пер. с англ. — М.: ИНФРА-М, 1999. — 402 с. — ISBN 8-86225-458-7.
• Кремер Н. Ш., Путко Б. А. Эконометрика. — М.: Юнити-Дана, 2003—2004. — 311 с. — ISBN 8-86225-458-7.
• Магнус Я. Р., Катышев П. К., Пересецкий А. А. Эконометрика. Начальный курс. — М.: Дело, 2007. — 504 с. — ISBN 978-5-7749-0473-0.
• Эконометрика. Учебник / Под ред. Елисеевой И. И. — 2-е изд. — М.: Финансы и статистика, 2006. — 576 с. — ISBN 5-279-02786-3.
• Александрова Н. В. История математических терминов, понятий, обозначений: словарь-справочник. — 3-е изд.. — М.: ЛКИ, 2008. — 248 с. — ISBN 978-5-382-00839-4.
• Витковский В. В. Наименьшие квадраты // Энциклопедический словарь Брокгауза и Ефрона : в 86 т. (82 т. и 4 доп.). — СПб., 1890—1907.
• Митин И. В., Русаков В. С. Анализ и обработка экспериментальных данных. — 5-е издание. — 24 с.

Ссылки[править | править код]

• Метод наименьших квадратов онлайн для зависимости y = a + bx с вычислением погрешностей коэффициентов и оцениванием автокорреляции.

Метод наименьших квадратов (МНК) — это статистическая процедура для довольно точного прогнозирования поведения зависимых переменных.

Например, можно понять, как будет меняться товарооборот (значение “y”) сети магазинов с изменением размеров торговой площади (значение “x”).

Суть МНК — из всех линейных функций найти наилучшее приближение к реальности. Это можно сделать путём поиска функции с наименьшим отклонением (точнее по процессу МНК: поиск минимальной суммы квадратов отклонений значений y (игрек) от полученного уравнения регрессии).

Решение МНК

Мы ищем уравнение линейной регрессии, которое выглядит так: y = ax + b

Где:

• y – зависимая переменная
• x – независимая переменная
• a – коэффициент (это также наклон/градиент линии)
• b – коэффициент (это также точка, где линия пересекает ось Y)

Метод 1

Шаги, которые мы будем делать для поиска y = ax + b (сейчас мы их все пройдём на примере):

Шаг 1: Для каждой точки (x, y) вычислить x² и xy.

Шаг 2: Суммировать все x, y, x² и xy, это даст нам Σx, Σy, Σx² и Σxy (если кто забыл, Σ означает “сумма”).

Шаг 3: Рассчитать наклон a по этой формуле:

, где N – количество данных

Шаг 4: Рассчитать значение числа b:

, где N – количество данных

Шаг 5: Подставить найденные числа по местам в уравнение (y = ax + b)

Пример

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

Размер (x)	Продажи (y)
2	4
3	6
5	8
7	10
9	12

Для наглядности, например, это магазины мороженого, и 2-метровая лавочка продаёт в месяц 4 тонны мороженого, 7 метровая — 10 тонн.

Шаг 1:

Сразу можно записать, что N = 5 (количество данных; т.е. всего у нас данные по пяти магазинам, ведь у нас 5 строк данных).

Для каждой точки (x, y) вычисляем x² и xy. Для этого, к уже существующим столбцам добавим ещё два: x² и xy.

• x² получим путём возведения x (Размер) в квадрат
• xy получим путём умножения одного на второе

x	y	x²	xy
2	4	2² = 4	2 × 4 = 8
3	6	3² = 9	3 × 6 = 18
5	8	25	40
7	10	49	70
9	12	81	108

Шаг 2: Суммировать все x, y, x² и xy, это даст нам Σx, Σy, Σx² и Σxy (складываем каждый столбик):

x	y	x²	xy
2	4	2² = 4	2 × 4 = 8
3	6	3² = 9	3 × 6 = 18
5	8	25	40
7	10	49	70
9	12	81	108
Σx = 26	Σy = 40	Σx² = 168	Σxy = 244

Шаг 3: Рассчитать a (наклон графика) по этой формуле:

, где N – количество данных

Помним, что N = 5, значит:

Шаг 4: Рассчитать значение числа b:

, где N – количество данных

Помним, что N = 5, значит:

Шаг 5: Подставить найденные числа по местам в уравнение

y = ax + b ⇒ y = 1,0976x + 2,29248

Готово!

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

Также мы можем использовать эту функцию, чтобы понять, как будут зависеть продажи фирмы от размера помещения. Например: руководство хочет открыть магазин размером в 11,5 м². Для этого подставляем 11,5 вместо x:

y = 1,0976x + 2,29248 ⇒ y = 1,0976 × 11,5 + 2,29248 = 14,91488

Ответ: этот магазин размером в 11,5 м² будет продавать около 15 тонн мороженого в месяц.

Метод 2

Мы продолжаем искать уравнение линейной регрессии, которое выглядит так: y = ax + b.

Используем тот же пример с сетью магазинов.

Размер (x)	Продажи (y)
2	4
3	6
5	8
7	10
9	12

Шаг 1: Опять суммируем все x, y, x² и xy, т.е. находим Σx, Σy, Σx² и Σxy (складываем каждый столбик):

x	y	x²	xy
2	4	2² = 4	2 × 4 = 8
3	6	3² = 9	3 × 6 = 18
5	8	25	40
7	10	49	70
9	12	81	108
Σx = 26	Σy = 40	Σx² = 168	Σxy = 244

Шаг 2: Записать вот такую систему уравнений (так мы будем искать параметры a и b):

Шаг 3: Помним, что N = 5. Таким образом, из нашего примера получаем систему:

Лучше конечно её переписать красиво:

Шаг 4: Решить систему.

Находим a = 1,0976; b = 2,29248; и ставим по местам в функцию (y = ax + b). Получается y = 1,0976x + 2,29248

Готово!

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

Узнайте также про Метод Крамера, Стандартное отклонение и Корреляции.