МНК: нормальная система уравнений и её решение

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

Что такое нормальная система МНК

Пусть у вас есть $n$ пар измерений $(x_i,\, y_i)$, и вы хотите аппроксимировать их линейной функцией $\hat{y} = a + bx$. Критерий МНК требует минимизировать сумму квадратов остатков:

$S(a,\,b) = \sum_{i=1}^{n} \bigl(y_i - a - b\,x_i\bigr)^2 \to \min.$

Необходимые условия минимума - равенство нулю частных производных по $a$ и $b$:

$\frac{\partial S}{\partial a} = -2\sum_{i=1}^{n}(y_i - a - b\,x_i) = 0,$

$\frac{\partial S}{\partial b} = -2\sum_{i=1}^{n}x_i(y_i - a - b\,x_i) = 0.$

Раскрывая суммы и перегруппировывая, приходим к нормальной системе МНК:

$\begin{cases} n\,a + \left(\sum x_i\right)b = \sum y_i, \\ \left(\sum x_i\right)a + \left(\sum x_i^2\right)b = \sum x_i y_i. \end{cases}$

Система содержит только пять статистик набора данных: $n$, $\sum x_i$, $\sum y_i$, $\sum x_i^2$, $\sum x_i y_i$ - всё, что нужно, чтобы перейти от точек к формулам.

Рассчитать для своей задачи

Как решить нормальную систему

Нормальная система для линейного случая - это $2 \times 2$. Решается подстановкой или формулами Крамера.

Обозначим суммы коротко:

$S_x = \sum_{i=1}^n x_i, \quad S_y = \sum_{i=1}^n y_i, \quad S_{x^2} = \sum_{i=1}^n x_i^2, \quad S_{xy} = \sum_{i=1}^n x_i y_i.$

Тогда определитель системы:

$D = n\,S_{x^2} - S_x^2.$

По формулам Крамера:

$b = \frac{n\,S_{xy} - S_x\,S_y}{D}, \qquad a = \frac{S_y - b\,S_x}{n}.$

Условие $D \neq 0$ выполняется всегда, если точки $x_i$ не все одинаковы (что очевидно для реального эксперимента).

Рассчитать для своей задачи

Геометрически каждое слагаемое $S(a,b)$ - это площадь квадрата со стороной, равной вертикальному отклонению точки от прямой. МНК выбирает прямую, при которой суммарная «площадь» всех таких квадратов наименьша.

Числовой пример: шаг за шагом

Возьмём 5 точек: $(1;\;2{,}1)$, $(2;\;3{,}9)$, $(3;\;6{,}2)$, $(4;\;8{,}0)$, $(5;\;9{,}8)$.

Шаг 1. Вычислим суммы:

Шаг 2. Составим нормальную систему ($n=5$):

$\begin{cases} 5a + 15b = 30{,}0, \\ 15a + 55b = 109{,}5. \end{cases}$

Шаг 3. Найдём определитель и решим:

$D = 5 \cdot 55 - 15^2 = 275 - 225 = 50,$

$b = \frac{5 \cdot 109{,}5 - 15 \cdot 30{,}0}{50} = \frac{547{,}5 - 450}{50} = \frac{97{,}5}{50} = 1{,}95,$

$a = \frac{30{,}0 - 1{,}95 \cdot 15}{5} = \frac{30{,}0 - 29{,}25}{5} = \frac{0{,}75}{5} = 0{,}15.$

Итог: $\hat{y} = 0{,}15 + 1{,}95\,x$.

Коэффициент детерминации R^2

После нахождения коэффициентов естественно проверить, насколько хорошо прямая описывает данные. Это делает коэффициент детерминации $R^2$:

$R^2 = 1 - \frac{SS_\text{res}}{SS_\text{tot}}, \quad SS_\text{res} = \sum(y_i - \hat{y}_i)^2,\quad SS_\text{tot} = \sum(y_i - \bar{y})^2,$

где $\bar{y} = S_y/n$ - среднее значение $y$. Значение $R^2 \in [0;\,1]$: $R^2 = 1$ означает идеальное совпадение прямой с данными, $R^2 \approx 0$ - полное отсутствие линейной зависимости.

Для нашего примера: $\bar{y} = 6{,}0$,

$SS_\text{tot} = (2{,}1-6)^2 + \ldots + (9{,}8-6)^2 = 33{,}8, \quad SS_\text{res} \approx 0{,}125,$

$R^2 \approx 1 - \frac{0{,}125}{33{,}8} \approx 0{,}996.$

Высокое $R^2$ подтверждает: линейная модель описывает данные превосходно.

Откуда берётся прямая МНК: матричная интерпретация

Запишем задачу компактно. Пусть $\mathbf{y} = (y_1, \ldots, y_n)^T$ - вектор наблюдений, $\mathbf{c} = (a, b)^T$ - вектор неизвестных коэффициентов, а матрица плана:

$A = \begin{pmatrix} 1 & x_1 \\ 1 & x_2 \\ \vdots & \vdots \\ 1 & x_n \end{pmatrix}.$

Тогда модель - $A\mathbf{c} \approx \mathbf{y}$, и критерий МНК - $\|A\mathbf{c} - \mathbf{y}\|^2 \to \min$. Минимум достигается в точке, где $A\mathbf{c} - \mathbf{y}$ ортогонален столбцам $A$, то есть $A^T(A\mathbf{c} - \mathbf{y}) = \mathbf{0}$. Это и есть нормальная система в матричной форме:

$A^T A\,\mathbf{c} = A^T \mathbf{y}.$

Перемножьте $A^T A$ вручную - получите в точности матрицу нормальной системы из предыдущего раздела: элемент $(1,1) = n$, $(1,2) = S_x$, $(2,2) = S_{x^2}$ и т. д. Матричная запись удобна для обобщений: достаточно изменить $A$, чтобы перейти от линейной модели к любому другому набору базисных функций.

Нормальная система для полиномиальной модели

МНК обобщается на произвольный полином $\hat{y} = a_0 + a_1 x + a_2 x^2 + \ldots + a_k x^k$. Нормальная система тогда содержит $k+1$ уравнение и строится по тому же принципу - приравниваем нулю частные производные по каждому коэффициенту:

$\begin{pmatrix} n & S_x & S_{x^2} \\ S_x & S_{x^2} & S_{x^3} \\ S_{x^2} & S_{x^3} & S_{x^4} \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ a_2 \end{pmatrix} = \begin{pmatrix} S_y \\ S_{xy} \\ S_{x^2 y} \end{pmatrix}.$

Матрица нормальной системы называется матрицей Грама (или матрицей информации). Для практических расчётов с $k \geq 2$ удобнее использовать матричную форму $A^T A\,\mathbf{c} = A^T \mathbf{y}$, где $A$ - матрица плана. Именно в таком виде МНК реализован во всех пакетах численных расчётов.

Практические применения МНК

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

• Физический эксперимент. Измерение зависимости сопротивления от температуры $R(T) = R_0(1 + \alpha T)$: это линейная модель, два параметра $R_0$ и $\alpha$ находятся нормальной системой.
• Графическая поверка закона. Если линеаризованный вид закона Хука $(F = k\Delta l)$ дать на вход МНК, наклон прямой сразу даст жёсткость пружины $k$ с наименьшей погрешностью.
• Демография и экономика. Линейный тренд временного ряда (выручка, население) - классическое применение нормальной системы с $x_i$ = номер года.
• Химия. Калибровочная кривая спектрофотометра: оптическая плотность линейно зависит от концентрации, МНК восстанавливает эту зависимость из серии стандартных образцов.

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

Частые ошибки

• Забывают перевести $x_i$ в числа при таблично заданных данных (например, года или даты). Нормальная система работает только с числами.
• Делят суммы на $n$ до подстановки. Систему нужно составлять именно из сырых сумм $\sum x_i$, $\sum x_i^2$ и т. д., а не из средних - деление вносит ошибку в уравнения.
• Путают $SS_\text{res}$ и $SS_\text{tot}$ при расчёте $R^2$. $SS_\text{tot}$ - разброс относительно среднего $\bar{y}$, $SS_\text{res}$ - относительно прямой. Поменять местами - получить $R^2 > 1$ или отрицательное.
• Ставят прочерк при $D = 0$. Если все $x_i$ одинаковы, определитель обращается в ноль и система вырождена. В таком случае надо менять набор данных, а не искать другую формулу.
• Не проверяют остатки на систематику. Высокое $R^2$ не означает правильность модели: если остатки образуют дугу, линейная модель неправильна, нужна парабола.

FAQ

Почему нормальная система называется «нормальной»? Название исторически связано с понятием нормали в геометрии: условие минимума суммы квадратов эквивалентно тому, что вектор остатков $\mathbf{e} = \mathbf{y} - A\mathbf{c}$ ортогонален (нормален) столбцам матрицы плана $A$. Именно это ортогональное условие и записывается в виде системы $A^T A\,\mathbf{c} = A^T \mathbf{y}$.

Можно ли применять МНК, если ошибки $y_i$ неодинаковы? Да, это взвешенный МНК: каждое наблюдение умножается на вес $w_i = 1/\sigma_i^2$. Нормальная система приобретает вид $A^T W A\,\mathbf{c} = A^T W \mathbf{y}$, где $W = \mathrm{diag}(w_i)$. При равных весах взвешенный МНК совпадает с обычным.

Что делать, если данных мало (например, $n = 3$) и нужна парабола? При числе параметров модели равном числу точек система имеет единственное решение, и $R^2 = 1$ автоматически - прямая (или парабола) проходит через все точки. Это не МНК-оптимизация, а точная интерполяция. Смысл МНК как аппроксимации появляется только при $n > k+1$.

Коротко

Нормальная система МНК получается из условия минимума суммы квадратов остатков: приравнивая нулю частные производные по всем коэффициентам, мы получаем линейную систему с матрицей Грама. Для линейной модели $\hat{y} = a + bx$ система имеет размер $2 \times 2$ и решается явными формулами через пять статистик набора данных. Качество подгонки оценивается коэффициентом детерминации $R^2$: чем он ближе к 1, тем лучше модель описывает разброс данных.