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

Метод найменших квадратів — метод знаходження наближеного розв'язку надлишково-визначеної системи. Часто застосовується в регресійному аналізі. На практиці найчастіше використовується лінійний метод найменших квадратів, що використовується у випадку системи лінійних рівнянь. Зокрема важливим застосуванням у цьому випадку є оцінка параметрів у лінійній регресії, що широко застосовується у математичній статистиці і економетриці.

Мотиваційний приклад

Графік точок даних (червоним), лінія найменших квадратів (синім) і відстані (зеленим)

Нехай в результаті деякого досліду отримано чотири ${\displaystyle (x,y)}$ точки даних: ${\displaystyle (1,6),}$ ${\displaystyle (2,5),}$ ${\displaystyle (3,7)}$ і ${\displaystyle (4,10)}$ (на малюнку ліворуч позначені червоним). Потрібно знайти пряму ${\displaystyle y=\beta _{1}+\beta _{2}x}$, яка найкраще підходить для цих точок. Інакше кажучи, ми хотіли б знайти числа ${\displaystyle \beta _{1}}$ і ${\displaystyle \beta _{2}}$, які приблизно розв'язують надвизначену лінійну систему

${\displaystyle {\begin{alignedat}{3}\beta _{1}+1\beta _{2}&&\;=\;&&6&\\\beta _{1}+2\beta _{2}&&\;=\;&&5&\\\beta _{1}+3\beta _{2}&&\;=\;&&7&\\\beta _{1}+4\beta _{2}&&\;=\;&&10&\\\end{alignedat}}}$

чотирьох рівнянь з двома невідомими в деякому найкращому сенсі.

Підхід найменших квадратів розв'язання цієї проблеми полягає у спробі зробити якомога меншою суму квадратів похибок між правою і лівою сторонами цієї системи, тобто необхідно знайти мінімум функції

${\displaystyle {\begin{aligned}S(\beta _{1},\beta _{2})=&\left[6-(\beta _{1}+1\beta _{2})\right]^{2}+\left[5-(\beta _{1}+2\beta _{2})\right]^{2}\\&+\left[7-(\beta _{1}+3\beta _{2})\right]^{2}+\left[10-(\beta _{1}+4\beta _{2})\right]^{2}.\end{aligned}}}$

Мінімум визначають через обчислення часткової похідної від ${\displaystyle S(\beta _{1},\beta _{2})}$ щодо ${\displaystyle \beta _{1}}$ і ${\displaystyle \beta _{2}}$ і прирівнюванням їх до нуля

${\displaystyle {\frac {\partial S}{\partial \beta _{1}}}=0=8\beta _{1}+20\beta _{2}-56}$

${\displaystyle {\frac {\partial S}{\partial \beta _{2}}}=0=20\beta _{1}+60\beta _{2}-154.}$

Це приводить нас до системи з двох рівнянь і двох невідомих, які називаються нормальними рівняннями. Роз'язком СЛАР будуть

${\displaystyle \beta _{1}=3.5}$

${\displaystyle \beta _{2}=1.4}$,

звідки отримуємо ${\displaystyle y=3.5+1.4x}$, що є рівнянням прямої, яка проходить найближче до поданих чотирьох точок. Мінімальна сума квадратів похибок є ${\displaystyle S(3.5,1.4)=1.1^{2}+(-1.3)^{2}+(-0.7)^{2}+0.9^{2}=4.2.}$

Результат підгонки сукупності спостережень ${\displaystyle (x_{i},y_{i})}$ (червоним) квадратичною функцією ${\displaystyle y=\beta _{1}+\beta _{2}x+\beta _{3}x^{2}\,}$ (синім). У лінійних найменших квадратах функція не повинна бути лінійною у своєму аргументі ${\displaystyle x,}$ а лише щодо своїх параметрів ${\displaystyle \beta _{j},}$ які треба визначити для отримання найкращого результату

Використання квадратичної моделі

Важливо, що у методі лінійних найменших квадратів ми не обмежені використанням прямої як моделі як у попередньому прикладі. Наприклад, ми могли вибрати обмежену квадратичну модель ${\displaystyle y=\beta _{1}x^{2}}$.^[1] Ця модель все ще лінійна в сенсі параметру ${\displaystyle \beta _{1}}$, отже ми все ще можемо здійснювати той самий аналіз, будуючи систему рівнянь з точок даних:

${\displaystyle {\begin{alignedat}{2}6&&\;=\beta _{1}(1)^{2}\\5&&\;=\beta _{1}(2)^{2}\\7&&\;=\beta _{1}(3)^{2}\\10&&\;=\beta _{1}(4)^{2}\\\end{alignedat}}}$

Часткові похідні щодо параметрів (цього разу лише одного) так само обчислюються і прирівнюються до 0:

${\displaystyle {\frac {\partial S}{\partial \beta _{1}}}=0=708\beta _{1}-498}$

Розв'язок отриманого рівняння:

${\displaystyle \beta _{1}=0.703,}$

що призводить до визначення найбільш підходящої моделі ${\displaystyle y=0.703x^{2}}$

Лінійний випадок

Одна незалежна змінна

Нехай маємо лінійну регресію зі скалярною змінною x:

${\displaystyle y=x\beta _{1}+\beta _{0},}$

а також вибірку початкових даних ${\displaystyle (y_{i},x_{i})}$ розміру M. Тоді

${\displaystyle \beta _{0}={\frac {1}{M}}\sum _{i}y_{i}-{\frac {\beta _{1}}{M}}\sum _{i}x_{i},\beta _{1}={\frac {M\sum _{i}x_{i}y_{i}-\sum _{i}x_{i}\sum _{i}y_{i}}{M\sum _{i}x_{i}^{2}-(\sum _{i}x_{i})^{2}}}}$

Множинна регресія (випадок багатьох незалежних змінних)

Для надлишково-визначеної системи m лінійних рівнянь з n невідомими ${\displaystyle \beta _{j},\quad (m>n):}$

${\displaystyle \sum _{j=1}^{n}X_{ij}\beta _{j}=y_{i},\quad i={\overline {1,m}},\quad j={\overline {1,n}}}$

чи в матричній формі запису:

${\displaystyle X{\boldsymbol {\beta }}=\mathbf {y} ,}$

зазвичай не існує точного розв'язку, і потрібно знайти такі β, які мінімізують наступну норму:

${\displaystyle {\underset {\boldsymbol {\beta }}{\operatorname {arg\,min} }}\,\sum _{i=1}^{m}\left|y_{i}-\sum _{j=1}^{n}X_{ij}\beta _{j}\right|^{2}={\underset {\boldsymbol {\beta }}{\operatorname {arg\,min} }}\,{\big \|}\mathbf {y} -X{\boldsymbol {\beta }}{\big \|}^{2}.}$

Такий розв'язок завжди існує і він є єдиним:

${\displaystyle {\hat {\boldsymbol {\beta }}}=(X^{\top }X)^{-1}X^{\top }\mathbf {y} }$

хоч дана формула не є ефективною через необхідність знаходити обернену матрицю.

Виведення формули

Значення ${\displaystyle S=\sum _{i=1}^{m}\left|y_{i}-\sum _{j=1}^{n}X_{ij}\beta _{j}\right|^{2}}$ досягає мінімуму в точці в якій похідна по кожному параметру рівна нулю. Обчислюючи ці похідні одержимо:

${\displaystyle {\frac {\partial S}{\partial \beta _{j}}}=2\sum _{i}r_{i}{\frac {\partial r_{i}}{\partial \beta _{j}}}=0\ (j=1,2,\dots ,n)}$

де використано позначення ${\displaystyle r_{i}=y_{i}-\sum _{j=1}^{n}X_{ij}\beta _{j}.}$

Також виконуються рівності:

${\displaystyle {\frac {\partial r_{i}}{\partial \beta _{j}}}=-X_{ij}.}$

Підставляючи вирази для залишків і їх похідних одержимо рівність:

${\displaystyle {\frac {\partial S}{\partial \beta _{j}}}=-2\sum _{i=1}^{m}X_{ij}\left(y_{i}-\sum _{k=1}^{n}X_{ik}\beta _{k}\right)=0.}$

Дану рівність можна звести до вигляду:

${\displaystyle \sum _{i=1}^{m}\sum _{k=1}^{n}X_{ij}X_{ik}{\hat {\beta }}_{k}=\sum _{i=1}^{m}X_{ij}y_{i}\ (j=1,2,\dots ,n)\,}$

або в матричній формі:

${\displaystyle (\mathbf {X} ^{\top }\mathbf {X} ){\hat {\boldsymbol {\beta }}}=\mathbf {X} ^{\top }\mathbf {y} .}$

Числові методи для обчислення розв'язку

Якщо матриця ${\displaystyle \ X^{\top }X}$ є невиродженою та додатноозначеною, тобто має повний ранг, тоді система може бути розв'язана за допомогою розкладу Холецького ${\displaystyle X^{\top }X=R^{\top }R}$, де ${\displaystyle R}$ — верхня трикутна матриця.

${\displaystyle R^{\top }R{\hat {\boldsymbol {\beta }}}=X^{\top }\mathbf {y} .}$

Розв'язок отримаємо в два кроки:

1. Отримаємо ${\displaystyle \mathbf {z} }$ з рівняння ${\displaystyle R^{\top }\mathbf {z} =X^{\top }\mathbf {y} ,}$
2. Підставимо і отримаємо ${\displaystyle {\hat {\boldsymbol {\beta }}}}$ з ${\displaystyle R{\hat {\boldsymbol {\beta }}}=\mathbf {z} .}$

В обох випадках використовуються властивості трикутної матриці.

Статистичні властивості

Одним із найважливіших застосувань лінійного МНК є оцінка параметрів лінійної регресії. Для заданого набору даних ${\displaystyle \{y_{i},\,x_{i1},\ldots ,x_{ip}\}_{i=1}^{n}}$ будується модель:

${\displaystyle y_{i}=\beta _{1}x_{i1}+\cdots +\beta _{p}x_{ip}+\varepsilon _{i}=x'_{i}\beta +\varepsilon _{i},\qquad i=1,\ldots ,n,}$

або в матричній формі:

${\displaystyle y=X\beta +\varepsilon ,\,}$

де:

${\displaystyle y={\begin{pmatrix}y_{1}\\y_{2}\\\vdots \\y_{n}\end{pmatrix}},\quad X={\begin{pmatrix}x'_{1}\\x'_{2}\\\vdots \\x'_{n}\end{pmatrix}}={\begin{pmatrix}x_{11}&\cdots &x_{1p}\\x_{21}&\cdots &x_{2p}\\\vdots &\ddots &\vdots \\x_{n1}&\cdots &x_{np}\end{pmatrix}},\quad \beta ={\begin{pmatrix}\beta _{1}\\\vdots \\\beta _{p}\end{pmatrix}},\quad \varepsilon ={\begin{pmatrix}\varepsilon _{1}\\\varepsilon _{2}\\\vdots \\\varepsilon _{n}\end{pmatrix}}.}$

В цих формулах ${\displaystyle \beta }$ — вектор параметрів, які оцінюються, наприклад, за допомогою методу найменших квадратів, а ${\displaystyle \varepsilon }$ — вектор випадкових змінних.

У класичній моделі множинної лінійної регресії приймаються такі умови:

• ${\displaystyle y_{i}=\beta _{0}\beta _{1}x_{i1}+\cdots +\beta _{p}x_{ip}+\varepsilon _{i}=x'_{i}\beta +\varepsilon _{i},\qquad i=1,\ldots ,n,}$
• ${\displaystyle \operatorname {E} [\,\varepsilon _{i}]=0.}$
• ${\displaystyle \operatorname {E} [\,\varepsilon _{i}\varepsilon _{j}]={\begin{cases}\sigma ^{2}&i=j\\0&i\neq j\end{cases}}}$

тобто випадкові змінні є гомоскедастичними і між ними відсутня будь-яка залежність.

• Ранг матриці X рівний p + 1, тобто між пояснюючими змінними відсутня лінійна залежність.

Для такої моделі оцінка ${\displaystyle {\hat {\boldsymbol {\beta }}}}$ одержана методом найменших квадратів володіє властивостями:

• Незміщеність. Оцінка ${\displaystyle {\hat {\boldsymbol {\beta }}}}$ є незміщеною, тобто ${\displaystyle \operatorname {E} [\,{\hat {\beta }}\,|X\,]=\beta .}$ Справді:

${\displaystyle \operatorname {E} [\,{\hat {\beta }}]=\operatorname {E} {\Big [}(X^{\top }X)^{-1}X^{\top }(X\beta +\varepsilon ){\Big ]}=\beta +\operatorname {E} {\Big [}(X^{\top }X)^{-1}X^{\top }\varepsilon {\Big ]}=\beta +(X^{\top }X)^{-1}X^{\top }\operatorname {E} (\varepsilon )=\beta }$

• Коваріаційна матриця оцінки ${\displaystyle {\hat {\boldsymbol {\beta }}}}$ рівна:

${\displaystyle \operatorname {Var} [\,{\hat {\beta }}\,]=\sigma ^{2}(X^{\top }X)^{-1}.}$

Це випливає з того, що ${\displaystyle \operatorname {Var} [\,Y\,]=\operatorname {Var} [\,\varepsilon \,]}$ і

${\displaystyle \operatorname {E} [\,{\hat {\beta }}]=\operatorname {Var} [\,(X^{\top }X)^{-1}X^{\top }Y\,]=(X^{\top }X)^{-1}X^{\top }\operatorname {Var} [\,Y\,]X(X^{\top }X)^{-1}=}$

${\displaystyle =\sigma ^{2}(X^{\top }X)^{-1}(X^{\top }X)(X^{\top }X)^{-1}=\sigma ^{2}(X^{\top }X)^{-1}}$

• Ефективність. Згідно з теоремою Гауса — Маркова оцінка, що одержана МНК, є найкращою лінійною незміщеною оцінкою.
• Змістовність. При доволі слабких обмеженнях на матрицю X метод найменших квадратів є змістовним, тобто при збільшенні розміру вибірки, оцінка за імовірністю прямує до точного значення параметру. Однією з достатніх умов є наприклад прямування найменшого власного значення матриці ${\displaystyle (X^{\top }X)}$ до безмежності при збільшенні розміру вибірки.
• Якщо додатково припустити нормальність змінних ${\displaystyle \varepsilon ,}$ то оцінка МНК має розподіл:

${\displaystyle {\hat {\beta }}\ \sim \ {\mathcal {N}}{\big (}\beta ,\ \sigma ^{2}(X^{\top }X)^{-1}{\big )}}$

В математичному моделюванні

Нехай ми маємо вибірку початкових даних ${\displaystyle f(x_{i})=y_{i}\ i={\overline {1..n}}}$. Функція ${\displaystyle f}$ — невідома.

Якщо ми знаємо приблизний вигляд функції ${\displaystyle f(x)}$, то задамо її у вигляді функціоналу ${\displaystyle F(x_{i},a_{0},\ldots ,a_{m})\approx y_{i}}$, де ${\displaystyle a_{0},\ldots ,a_{m}}$ — невідомі константи.

Нам потрібно мінімізувати відмінності між ${\displaystyle F}$ та ${\displaystyle f}$. Для цього беруть за міру суму квадратів різниць значень цих функцій у всіх точках ${\displaystyle x_{i}}$ і її мінімізують (тому метод так і називається):

${\displaystyle I(a_{0},\ldots ,a_{m})=\sum _{i=0}^{n}(y_{i}-F(x_{i},a_{0},\ldots ,a_{m}))^{2}\to \min }$

Коефіцієнти ${\displaystyle a_{j}}$ в яких така міра мінімальна знаходять з системи:

${\displaystyle {\begin{cases}\displaystyle {\frac {\partial I(a_{0},\ldots ,a_{m})}{\partial a_{0}}}=0\\\ldots \\\displaystyle {\frac {\partial I(a_{0},\ldots ,a_{m})}{\partial a_{m}}}=0\end{cases}}}$