W poprzednim artykule pisałem właśnie o regresji liniowej i jej zastosowaniu. Niniejszy artykuł jest kontynuacją tej pracy, w którym formalnie wyprowadzimy estymatory metody najmniejszych kwadratów parametrów modelu regresji liniowej.
Załóżmy, że dysponujemy \(n\) obserwacjami \(y_1,y_2,\dots,y_n\) zmiennej \(y\) oraz odpowiadającymi im obserwacjami \(x_1,x_2,\dots,x_n\) zmiennej \(x\). Wtedy model zależności w regresji liniowej dla powyższych obserwacji jest postaci
\begin{equation}\label{model}
y_i=ax_i+b+e_i,\quad i=1,2,\dots,n,
\end{equation}
przy czym \(a\neq 0\), \(b\in R\) są nieznanymi parametrami, a \(e_i\) to błędy losowe. Oszacowania parametrów \(a\) i \(b\) uzyskane metodą najmniejszych kwadratów mają postać
\begin{equation}\label{estymatory}
\hat{a} =\frac{y_1(x_1-\bar{x})+y_2(x_2-\bar{x})+\dots+y_n(x_n-\bar{x})}{(x_1-\bar{x})^2+(x_2-\bar{x})^2+\dots+(x_n-\bar{x})^2},\quad
\hat{b}=\bar{y}-\bar{x}\hat{a},
\end{equation}
przy czym \(\bar{x}=\frac{1}{n}(x_1+x_2+\dots+x_n)\) i \(\bar{y}=\frac{1}{n}(y_1+y_2+\dots+y_n)\) są średnimi arytmetycznymi z obserwacji zmiennych \(x\) i \(y\), odpowiednio. Otrzymujemy zatem prostą postaci \(y=\hat{a}x+\hat{b}\) zwaną prostej regresji opisującą zależność między zmiennymi \(x\) i \(y\).
Wyprowadzenie estymatorów metody najmniejszych kwadratów parametrów \(a\) i \(b\) podanych we wzorach \eqref{estymatory} rozpoczniemy od przedstawienia modelu \eqref{model} w postaci wykorzystującej macierze. Poniżej przedstawimy kilka potrzebnych informacji o macierzach, lecz w trosce o zachowanie przejrzystości artykułu, ograniczamy formalizm i szczegółowe definicje do minimum.
Macierz to pewna funkcja, której jednak nie będziemy tutaj dokładnie definiować. Wystarczy nam intuicyjne zrozumienie macierzy jako tablicy liczb. Przykładowa macierz o czterech wierszach i trzech kolumnach jest następująca
\begin{equation}\label{przmacierzy}
\left[\begin{array}{ccc}
3&1&5\\
5&3&7\\
4&4&1\\
6&2&5\\
\end{array}\right].
\end{equation}
Liczby, które występują w macierzy nazywamy jej elementami. Elementami powyższej macierzy są liczby \(3,1,5,5,3,7,4,4,1,6,2,5\). Macierzy o tej samej liczbie wierszy i kolumn możemy dodawać i odejmować element po elemencie, np.
\[
\left[\begin{array}{ccc}
3&1&5\\
5&3&7\\
4&4&1\\
6&2&5\\
\end{array}\right]+\left[\begin{array}{ccc}
5&2&6\\
1&4&4\\
7&3&5\\
5&2&6\\
\end{array}\right]=
\left[\begin{array}{ccc}
3+5&1+2&5+6\\
5+1&3+4&7+4\\
4+7&4+3&1+5\\
6+5&2+2&5+6\\
\end{array}\right]
=
\left[\begin{array}{ccc}
8&3&11\\
6&7&11\\
11&7&6\\
11&4&11\\
\end{array}\right].
\]
Ponadto, macierze możemy również mnożyć. Aby wymnożyć dwie macierze przez siebie, mnożna musi mieć taką samą liczbę kolumn jak liczba wierszy mnożnika. Mnożenie macierzy nie jest mnożeniem element po elemencie jak w przypadku dodawania. Niech
\[
A=\left[ \begin{array}{cccc}
a_{11}&a_{12}&\cdots&a_{1n}\\
a_{21}&a_{22}&\cdots&a_{2n}\\
\vdots&\vdots&\ddots&\vdots\\
a_{m1}&a_{m2}&\cdots&a_{mn}
\end {array}\right],
\quad B=
\left[ \begin{array}{cccc}
b_{11}&b_{12}&\cdots&b_{1k}\\
b_{21}&b_{22}&\cdots&b_{2k}\\
\vdots&\vdots&\ddots&\vdots\\
b_{n1}&b_{n2}&\cdots&b_{nk}
\end {array}\right].
\]
Macierze \(A\) i \(B\) mnożymy według wzoru
\[
AB=\left[ \begin{array}{cccc}
c_{11}&c_{12}&\cdots&c_{1k}\\
c_{21}&c_{22}&\cdots&c_{2k}\\
\vdots&\vdots&\ddots&\vdots\\
c_{m1}&c_{m2}&\cdots&c_{mk}
\end {array}\right],
\]
przy czym
\[
c_{ij}=a_{i1}b_{1j}+a_{i2}b_{2j}+\dots+a_{in}b_{nj}.
\]
Będzie nam potrzebna jeszcze jedna operacja na macierzach, a mianowicie transpozycja macierzy. Macierz transponowana macierzy \(A\) (oznaczana \(A’\)) to taka, której wiersze są kolumnami macierzy \(A\). Mówiąc bardziej zrozumiale, zamieniamy tutaj wiersze na kolumny. Macierzą transponowaną macierzy \eqref{przmacierzy} jest macierz
\[
\left[\begin{array}{cccc}
3&5&4&6\\
1&3&4&2\\
5&7&1&5\\
\end{array}\right].
\]
Możemy teraz przejść do przedstawienia modelu \eqref{model} w postaci macierzowej. Zdefiniujmy wpierw macierze
\begin{equation}\label{yXbetae}
Y=\left[\begin{array}{c}
y_1\\
y_2\\
\vdots\\
y_n\\
\end{array}\right],\quad X=\left[\begin{array}{cc}
1&x_1\\
1&x_2\\
\vdots&\vdots\\
1&x_n\\
\end{array}\right],\quad\beta=\left[\begin{array}{c}
b\\
a\\
\end{array}\right],\quad E=\left[\begin{array}{c}
e_1\\
e_2\\
\vdots\\
e_n\\
\end{array}\right]
\end{equation}
oznaczające odpowiednio macierze obserwacji zmiennej \(y\), obserwacji zmiennej \(x\) wraz z kolumną samych jedynek odpowiadającą wyrazowi wolnemu \(b\), parametrów \(a\) i \(b\), oraz macierz błędów \(e_i,i=1,2,\dots,n\). Przy takich oznaczeniach model \eqref{model} możemy zapisać wzorem
\[
Y=X\beta+E.
\]
To przedstawienie modelu \eqref{model} pozwala nam skorzystać ze znanego wzoru na estymatory metody najmniejszych kwadratów parametrów będących elementami macierzy \(\beta\), przedstawiającego się następująco
\[
\hat{\beta}=(X’X)^{-1}X’Y,
\]
przy czym \((X’X)^{-1}\) oznacza odwrotność macierzy \(X’X\) (za chwilę powiemy o niej więcej). Podstawiając do tego wzoru macierze \(X\) i \(Y\) podane w \eqref{yXbetae}, wyznaczymy estymatory parametrów \(a\) i \(b\) zaprezentowane równaniami \eqref{estymatory}. Korzystając z podanych powyżej definicji, najpierw transponujemy macierz \(X\):
\[
X’=\left[\begin{array}{cccc}
1&1&\dots&1\\
x_1&x_2&\dots&x_n\\
\end{array}\right],
\]
a następnie wymnażamy przez siebie macierze \(X’\) i \(X\):1
\[
X’X=\left[\begin{array}{cccc}
1&1&\dots&1\\
x_1&x_2&\dots&x_n\\
\end{array}\right]\left[\begin{array}{cc}
1&x_1\\
1&x_2\\
\vdots&\vdots\\
1&x_n\\
\end{array}\right]=\left[\begin{array}{cc}
n&\sum_{i=1}^nx_i\\
\sum_{i=1}^nx_i&\sum_{i=1}^nx_i^2\\
\end{array}\right]=\left[\begin{array}{cc}
n&n\bar{x}\\
n\bar{x}&\sum_{i=1}^nx_i^2\\
\end{array}\right].
\]
Teraz musimy odwrócić macierz \(X’X\), czyli znaleźć do niej macierz odwrotną. Macierz odwrotna \(A^{-1}\) do macierzy \(A\) to taka macierz, dla której \(A^{-1}A=AA^{-1}=I\), przy czym
\[
I=\left[ \begin{array}{cccc}
1&0&\cdots&0\\
0&1&\cdots&0\\
\vdots&\vdots&\ddots&\vdots\\
0&0&\cdots&1
\end{array}\right]
\]
jest macierzą jednostkową. W świecie macierzy, macierz odwrotna jest czymś w rodzaju liczby odwrotnej (\(x\) — liczba niezerowa, \(1/x\) — liczba odwrotna do liczby \(x\)) a macierz jednostkowa — liczby jeden. Zazwyczaj szukanie macierzy odwrotnej jest czasochłonne. Na szczęście, nam jest potrzebna macierz odwrotna do macierzy o dwóch wierszach i kolumnach (patrz \(X’X\)), którą łatwo znaleźć korzystając ze wzoru
\[
\left[\begin{array}{cc}
a&b\\
c&d\\
\end{array}\right]^{-1}=\left[\begin{array}{cc}
\frac{d}{ad-cb}&\frac{-b}{ad-cb}\\
\frac{-c}{ad-cb}&\frac{a}{ad-cb}\\
\end{array}\right].
\]
Mianownik ,,\(ad-cb\)” dla macierzy \(X’X\) liczymy następująco
\[
n\cdot \sum_{i=1}^nx_i^2-(n\bar{x})^2=n\left(\sum_{i=1}^nx_i^2-n\bar{x}^2\right)=n\sum_{i=1}^n(x_i-\bar{x})^2,
\]
przy czym ostatnią równość otrzymujemy w następujący sposób korzystając ze wzoru skróconego mnożenia
\begin{align*}
\sum_{i=1}^n(x_i-\bar{x})^2&=\sum_{i=1}^n(x_i^2-2x_i\bar{x}+\bar{x}^2)\\
&=\sum_{i=1}^nx_i^2-\sum_{i=1}^n2x_i\bar{x}+\sum_{i=1}^n\bar{x}^2=\sum_{i=1}^nx_i^2-2\bar{x}\sum_{i=1}^nx_i+\bar{x}^2\sum_{i=1}^n1\\
&=\sum_{i=1}^nx_i^2-2\bar{x}(n\bar{x})+n\bar{x}^2=\sum_{i=1}^nx_i^2-2n\bar{x}^2+n\bar{x}^2=\sum_{i=1}^nx_i^2-n\bar{x}^2.
\end{align*}
Zatem
\[
(X’X)^{-1}=\left[\begin{array}{cc}
\frac{\sum_{i=1}^nx_i^2}{n\sum_{i=1}^n(x_i-\bar{x})^2}&\frac{-n\bar{x}}{n\sum_{i=1}^n(x_i-\bar{x})^2}\\
\frac{-n\bar{x}}{n\sum_{i=1}^n(x_i-\bar{x})^2}&\frac{n}{n\sum_{i=1}^n(x_i-\bar{x})^2}\\
\end{array}\right].
\]
Teraz wyznaczamy \(X’Y\):
\[
X’Y=\left[\begin{array}{cccc}
1&1&\dots&1\\
x_1&x_2&\dots&x_n\\
\end{array}\right]\left[\begin{array}{c}
y_1\\
y_2\\
\vdots\\
y_n\\
\end{array}\right]=\left[\begin{array}{c}
\sum_{i=1}^ny_i\\
\sum_{i=1}^nx_iy_i\\
\end{array}\right]=\left[\begin{array}{c}
n\bar{y}\\
\sum_{i=1}^nx_iy_i\\
\end{array}\right].
\]
Podsumowując
\[
\hat{\beta}=\left[\begin{array}{c}
\hat{b}\\
\hat{a}\\
\end{array}\right]=\left[\begin{array}{cc}
\frac{\sum_{i=1}^nx_i^2}{n\sum_{i=1}^n(x_i-\bar{x})^2}&\frac{-n\bar{x}}{n\sum_{i=1}^n(x_i-\bar{x})^2}\\
\frac{-n\bar{x}}{n\sum_{i=1}^n(x_i-\bar{x})^2}&\frac{n}{n\sum_{i=1}^n(x_i-\bar{x})^2}\\
\end{array}\right]\left[\begin{array}{c}
n\bar{y}\\
\sum_{i=1}^nx_iy_i\\
\end{array}\right]=\left[\begin{array}{cc}
\frac{(\sum_{i=1}^nx_i^2)n\bar{y}-n\bar{x}\sum_{i=1}^nx_iy_i}{n\sum_{i=1}^n(x_i-\bar{x})^2}\\
\frac{-n^2\bar{x}\bar{y}+n\sum_{i=1}^nx_iy_i}{n\sum_{i=1}^n(x_i-\bar{x})^2}\\
\end{array}\right].
\]
Po przekształceniach, które zostawiamy Czytelnikowi jako ćwiczenie, otrzymujemy estymatory \(\hat{a}\) i \(\hat{b}\) podane w \eqref{estymatory}.
Przypisy
- Korzystamy ze skróconego zapisu sumowania za pomocą znaku sumy \(\sum\). Przykładowo, sumę \(a_1+a_2+\dots+a_n\) możemy w skrócie zapisać następująco \(\sum_{i=1}^na_i\).
Artykuł został sfinansowany dzięki wsparciu pozyskanemu przez Poznańską Fundację Matematyczną od Miasta Poznań na realizację projektu ,,Potęga matematyki''.