Metoda najmniejszych kwadratów

Metoda najmniejszych kwadratów

Metoda najmniejszych kwadratów, a właściwie metoda najmniejszej sumy kwadratów odchyleń, to jedna z najprostszych i najstarsza metoda rozwiązania problemu regresji. W problemie regresji mamy do czynienia z pomiarem dwóch (lub więcej) zmiennych, dokonywanym w sposób sparowa- ny, tzn. każdej wartości zmiennej przypisujemy jedną wartość zmiennej niezależnej i poszukujemy związku pomiędzy zmienną niezależną a zmienną zależną. Wybór zmiennej zależnej nie zawsze musi być jednoznaczny, ale w zwykłej metodzie najmniejszych kwadratów przyjmujemy, że zmienna niezależna jest wyznaczona bez żadnych odchyłek od wartości rzeczywistej, natomiast zmienna zależna takie odchyłki zawiera. Odchyłki, zwane też często (niewłaściwie) błędami, mogą być spo- wodowane niedokładnością pomiaru, wynikającą z właściwości miernika, albo wpływem czynników zakłócających pomiar, których nie możemy kontrolować.

W praktyce rozwiązanie problemu regresji polega na dopasowaniu parametrów funkcji do danych pomiarowych. Zabieg ten może służyć dwom celom:

I. Chcemy odkryć rodzaj zależności (chcemy dobrać wzór empiryczny, który najlepiej opisuje zmie- rzoną zależność).

II. Chcemy poznać parametry danej teoretycznej zależności poprzez dopasowanie ich do obserwo- wanych wartości.

Mając zbiór punktów uzyskanych wskutek obserwacji dla których zależność wielkości y od x opisuje funkcja y= f (x , p_1,p_2,…, p_k), której kształt (np. nachylenie, liczba punktów przegięcia, zakrzy- wienie, itp.) zależy od pewnych parametrów: p1, p2, … , pk.

0 1 2 3 4 5 6 7 8

0,0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9

Metoda najmniejszych kwadratów polega na tym, żeby tak dobrać te parametry, aby suma odleg- łości punktów pomiarowych od krzywej przebiegu funkcji była jak najmniejsza. Jeżeli krzywa ma skomplikowany przebieg, to wyznaczenie wzorów na odległość punktu od niej (odcinek prostopadły do krzywej) staje się trudne i dlatego w najprostszym podejściu poszukuje się minimum sumy odległości tylko w jednym wymiarze zmiennej zależnej i taka metoda nazywana jest zwyczajną metodą najmniejszych kwadratów.

W metodzie tej dla danego pomiaru i oblicza się różnicę pomiędzy teoretyczną wartością funkcji w punkcie xi, równą f(xi), a obserwowaną wartością yi. Na obrazku odległości te przedstawiają czar- ne kreski. Celem metody najmniejszych kwadratów jest wyznaczenie takiego zestawu parametrów funkcji, aby suma kwadratów tych różnic była minimalna.

0 1 2 3 4 5 6 7 8

0,0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9

Suma kwadratów odchyleń:

S = ∑

i

( y

− f ( x

, p

p

… , p

))

jest funkcją celu i można formalnie traktować ją jako funkcję parametrów p1, p2, … , pk, gdyż obserwowane wartości xi i yi, są zadane przez obserwację, a jej wartość zależy jedynie od wartości parametrów.

Wyznaczenie minimum funkcji sprowadza się więc do obliczenia dla każdego parametru pj różniczki cząstkowej ∂S/∂pj i przyrównanie jej do zera, w rezultacie czego dostajemy układ k równań na k niewiadomych (p1, p2, … , pk), co jednoznacznie wyznacza ich wartości. Nie zawsze jednak układ takich równań da się rozwiązać metodami analitycznymi, dlatego też w praktyce stosuje się często metody numeryczne.

W przypadku Excela mamy bezpośrednią możliwość wyznaczenia metodami numerycznymi para- metrów funkcji za pomocą linii trendu na wykresie dla 5 rodzajów zależności:

1) liniowej:

y=ax+b

2) logarytmicznej: y=a ln x+b; 3) potęgowej:

y=ax

4) wykładniczej (eksponencjalnej): y=a e^bx i

5) wielomianowej:

y=a

x

+a

x

+…+a

x+a

Wyznaczanie parametrów linii trendu

Wyznaczanie parametrów linii trendu zostało pokazane dla doboru wielomianu 5 stopnia w pliku mnk.xlsx, w zakładce mnk-wiel dokumentu. Dane to zmierzone widmo linii promieniowania, gdzie zmienna lambda to długość fali promieniowania wyrażona w Angstremach (1Å=10⁻¹⁰m), a Intesity to natężenie promieniowania w tej długości wyrażona w jednostkach umownych. Kolej- ność czynności jest następująca.

1. Tworzymy wykres zależności Y od X, w tym wypadku zależności Intesity od lambda.

2. Klikamy prawym klawiszem na serię danych, dla której chcemy wyznaczyć parametry linii tren- du (w tym wypadku mamy tylko jedną serię na wykresie), i z menu podręcznego wybieramy opcję Dodaj linię trendu.

3. W pojawiającym się oknie wybieramy rodzaj dobieranej funkcji. W przypadku wyboru typu wielomianowego, należy jeszcze wybrać stopień wielomianu (nie wiem, dlaczego w tej wersji Excela rubryka ta nosi nazwę Kolejność; w innych wersjach może to być jakaś inna nazwa, ale zawsze leży przy wyborze opcji wielomianowej). Dobrze jest jeszcze wybrać opcje: Wyświetl równanie na wykresie oraz Wyświetl wartości R-kwadrat na wykresie.

Wartość R² to współczynnik determinacji, który pokazuje jaka część zmienności zmiennej zależnej jest wyjaśniana przez regresję. Wyznaczone parametry (współczynniki) wielomianu możemy odczy- tać z wyświetlonego na wykresie równania. Wyświetloną na wykresie linię trendu można formato- wać, klikając na nią prawym klawiszem i wybierając z menu podręcznego odpowiednie elementy, które chcemy zmienić (ale są to oczywiście jedynie zabiegi kosmetyczne). Ja zmieniłem kolor na czerwony i zwiększyłem grubość linii.

Mimo bardzo dobrej wartości współczynnika R², wykres linii trendu pokazuje, że wielomian nie jest dobrym przybliżeniem tej zależności (proponuję sprawdzić inne dostępne linie trendu i inne stopnie wielomianu). Dobrze więc, gdy skądinąd możemy dobrać właściwszy przebieg funkcji. Możemy go ustalić na podstawie kształtu ogólnego przebiegu zależności, albo na podstawie teorii opisującej daną zależność, jeśli ją znamy.

Wyznaczenie parametrów funkcji własnej

W przypadku funkcji innej niż proponowane przez Excel możemy posłużyć się albo metodą wyrów- nywania, albo bezpośrednio metodą najmniejszych kwadratów.

Metoda wyrównywania polega na tym, że spośród funkcji o najbardziej podobnym przebiegu do za- leżności zmierzonej, takich że możliwe jest przekształcenie zmierzonych wielkości x i y w wielkości x′=φ(x,y) i y′=ψ(x,y), tak aby zależność pomiędzy y′ a x′ była liniowa, wybiera się kilka, a następ- nie sprawdza się, dla której zależność pomiędzy y′ a x′ jest najbliższa liniowej.

Jeżeli mamy jednak jakieś podstawy teoretyczne do wyboru postaci funkcji opisującej zależność, to lepiej posłużyć się bezpośrednio metodą najmniejszych kwadratów.

Takie właśnie wyznaczanie parametrów funkcji zadanej zostało pokazane w pliku mnk.xlsx, w za- kładce mnk-wiel dokumentu. Dane to, to samo, co poprzednio, zmierzone widmo linii promieniowa- nia, gdzie zmienna lambda to długość fali promieniowania wyrażona w Angstremach (1Å=10⁻¹⁰m), a Intesity to natężenie promieniowania w tej długości wyrażona w jednostkach umownych. Kształt profilu linii widmowej promieniowania atomu przy przejściu z jednego poziomu energetycznego na inny jest opisywany przez funkcję Lorenza:

I = A

1+ ( ^λ

^{Δ −λ 2} )

^+I

gdzie I oraz λ to nasze zmienne Intesity oraz lambda, natomiast parametr A to intensywność linii, parametr λ0 to położenie linii (środka jej profilu), Δ to szerokość połówkowa linii, a I0 to poziom tła.

Dokładne wartości parametrów A, λ0, Δ oraz I0 zostaną dobrane metodą najmniejszych kwadratów za pomocą minimalizacji funkcji celu (sumy kwadratów odchyleń) metodą numeryczną przy użyciu narzędzia Solver.

Kolejność czynności jest następująca.

1. Jak poprzednio, tworzymy wykres zależności Y od X, w tym wypadku zależności Intesity od lambda.

2. W osobnych komórkach wpisujemy początkowe wartości parametrów funkcji: A w komórce G1, λ0 w komórce G2, Δ w komórce G3, I0 w komórce G4.

Tutaj bardzo ważna uwaga! Jeżeli funkcja jest liniowa względem parametrów, tzn. postaci:

f (x )= p₁ f₁(x )+ p₂ f₂(x)+…+ p_k f_k(x)

to nie ma większego znaczenia, od jakich wartości parametrów startujemy, ale w ogólności nie zawsze możemy trafić na rozwiązanie, gdyż metody numeryczne poszukują minimum funkcji metodą kolejnych przybliżeń. Poszukiwanie takie może zakończyć się np. w miejscu jakiegoś mi- nimum lokalnego funkcji celu, albo nie doprowadzić do żadnego rozwiązania. Ważny jest więc odpowiedni dobór wartości startowych parametrów. Można je dobrać albo na podstawie jakichś rozważań teoretycznych, albo, jak w tym przypadku, na podstawie wykresu. Wybrane na tej podstawie wartości początkowe były następujące: A=700000, λ0=1052, Δ =1,5, I0=50000.

3. W kolumnie następnej (C) liczymy wartość funkcji, posługując się adresami bezwzględnymi ko- mórek z wartościami parametrów. Uwaga! Nie podajemy wartości parametrów we wzorze za pomocą liczb, gdyż traktujemy je w metodzie jako zmienne.

4. Nie jest to konieczne, ale warto po obliczeniu wartości funkcji, wprowadzić jej wartości na wy- kres jako dodatkową serię, gdyż będzie to pomocne przy ocenie, czy dobór wartości parame- trów jest właściwy. Najlepiej serię tę wykreślić w postaci linii. Rezultat widoczny poniżej.

5. W kolejnej kolumnie (D) liczymy kwadraty różnic pomiędzy wartością obserwowaną zmiennej zależnej (Intensity) a wartością funkcji.

6. W komórce G5 liczymy sumę tych kwadratów (sumę wartości z kolumny D). Celem metody jest zminimalizowanie jej wartości poprzez dobranie odpowiednich wartości parametrów funkcji.

7. Aby zminimalizować wartość sumy kwadratów odchyleń, należy użyć narzędzia Solver, które również znajduje się w zakładce Dane, w grupie Analiza. Ponieważ celem jest minimum sumy kwadratów odchyleń, w polu Ustaw cel: podajemy adres komórki (G5), gdzie obliczana jest ta suma, a następnie w polu Na: wybieramy opcję Min. W polu o nazwie Przez zmienianie komórek zmiennych: podajemy zakres komórek, w których podane są wartości parametrów funkcji (G1:G4), ponieważ jedynie poprzez ich zmianę mamy wpływ na wartość sumy kwadra- tów odchyleń oraz ponieważ to te parametry chcemy dopasować. Poniżej, w polu o nazwie Podlegających ograniczeniom:, możemy podać dodatkowe warunki ograniczające wartoś- ci parametrów. Źródłem tych dodatkowych warunków mogą być na przykład jakieś właściwości parametrów wynikające z teorii opisującej badaną zależność (np. w naszym przypadku inten- sywność linii A nie może być ujemna, tak samo długość fali czy szerokość profilu). W tym wypadku zostało to zaznaczone przez globalne wybranie opcji Ustaw wartości nieujemne dla zmiennych bez ograniczeń poniżej tego pola, gdyż w tym wypadku to ograniczenie dotyczy wszystkich parametrów. W polu poniżej o nazwie Wybierz metodę rozwiązania:

została wybrana metoda odpowiednia dla funkcji nieliniowych (Nieliniowa GRG).

8. Po naciśnięciu przycisku Rozwiąż pojawia się rozwiązanie wraz z oknem opisującym efekty kolejnych przybliżeń. Jak widać, chociażby z przebiegu linii funkcji na wykresie, rezultat jest cał- kowicie niepoprawny. Należy więc wybrać opcję Przywróć wartości pierwotne i wcisnąć przycisk OK.

9. W przypadku braku zadowalającego rozwiązania można przyjąć dwie strategie:

1) uściślić warunki ograniczające dla parametrów, jeśli mamy na ten temat jakieś wyobrażenie;

2) zmienić wartości startowe parametrów (niekoniecznie wszystkich).

W tym wypadku widać, że jedynym słabo określonym parametrem był parametr A, gdyż pozo- stałe były w miarę dobrze określone na podstawie obserwacji wykresu. Podejrzenie to potwier- dza również wynik pierwszego dopasowania, gdzie, jak widać, wartość tego parametru poszybo- wała do ogromnych, nierealistycznych wartości.

10. W drugim podejściu ustawiłem więc wartość początkową tego parametru na 500000, a pozosta- łe otrzymały te same wartości początkowe, co poprzednio. W wyniku uruchomieniu Solvera, tym razem rezultat był zadowalający, co widać w załączonym pliku.

11. Ponieważ liczba punktów jest na tyle mała, że wykres funkcji, zwłaszcza w części środkowej, jest zbyt łamany, wprowadziłem dodatkową serię na wykres, opartą na tabelce funkcji dla wartości lambda w zakresie od 1050 do 1055 co 0,01 (kolumny U i V). Dzięki temu otrzymałem gładki wykres funkcji. Nie jest to oczywiście konieczne i nic nie wnosi do rozwiązania, poza efektem estetycznym.

Jak to zostało wspomniane w przypadku zastosowania linii trendu na wykresie, jednym ze wskaźni- ków pozwalających określić jakość dopasowania parametrów jest współczynnik determinacji R². W przypadku posłużenia się linią trendu, jest on obliczany automatycznie, tutaj obliczymy go samo- dzielnie, korzystając z jego definicji. Współczynnik determinacji R² jest proporcją modelowej sumy kwadratów odchyleń SSM do całkowitej sumy kwadratów odchyleń SST, czyli R²=SSM/SST. SST to suma kwadratów różnic pomiędzy wartościami obserwowanymi zmiennej zależnej a ich wartością średnią. SSM to suma kwadratów różnic pomiędzy wartościami modelowanymi przez funkcję a ich wartością średnią zmiennej zależnej.

Obliczenia te są przeprowadzone na arkuszu R2.

1. Na arkusz skopiowałem wartości obserwowane Y, czyli Intensity do kolumny A, oraz wartości przewidywane przez model, czyli wartości funkcji (funkcja do kolumny B) już po ustaleniu war- tości parametrów za pomocą metody najmniejszych kwadratów. Żeby nie kopiować całości ar- kusza, a tylko wartości funkcji, przekopiowałem je jako same wartości liczbowe (pokazywałem jak to robić w poprzednich tematach).

2. Ponieważ chcę docelowo liczyć kwadraty odchyleń od wartości średniej obserwowanych Y, więc w komórce F1 obliczyłem za pomocą funkcji ŚREDNIA() wartość średnią liczb z kolumny A.

3. W kolumnie C (ST) obliczyłem kwadraty odchyleń obserwowanych wartości Y (Intensity) od wartości średniej, a w kolumnie D (SM) obliczyłem kwadraty odchyleń wartości modelowanych (funkcja) od wartości średniej.

4. W komórce F2 (SSM) obliczyłem całkowitą sumę kwadratów odchyleń (sumę z kolumny D), a w komórce F3 (SST) obliczyłem modelową sumę kwadratów odchyleń (sumę z kolumny C).

5. Ostatecznie w komórce F4 obliczyłem wartość współczynnika determinacji R² zgodnie z jego definicją. Jak widać, jest on całkiem duży i większy niż dla linii trendu wielomianowej.

Zadania do wykonania

1. Spróbować dopasować inne linie trendu i przeanalizować wyniki (wykres, otrzymane równanie i wartość R²).

2. Dla linii trendu wielomianowego sprawdzić, jak zmienia się wartość R² dla różnych stopni wielo- mianu.

3. Metodą najmniejszych kwadratów dopasować do obserwowanych wartości funkcję Gaussa. Jest to profil linii widma promieniowania, gdy dominujący wpływ na jej kształt mają efekty Dopplera.

Przyjąć funkcję w postaci:

I =A exp ( ^{− ( λ}

^{Δ −λ 2 )}

) ^{+ I}

gdzie I oraz λ to zmienne Intesity oraz lambda, parametr A to intensywność linii, parametr λ0

to położenie linii (środka jej profilu), Δ to szerokość połówkowa linii, a I0 to poziom tła.

Uwaga! Ta wersja Excela, w której pracowałem, ma tę nieprzyjemną wartość, że −(coś)^2 jest liczone jak (−(coś))^2, co jest błędem. Aby się przed tym zabezpieczyć, zamiast −(coś)^2 należy pisać −((coś)^2), co da właściwy wynik. Nie daję głowy, że w Państwa wersji tak będzie,

więc najlepiej sprawdzić to we własnej wersji Excela doświadczalnie. Przypominam, że funkcja exp(x) to e^x.

4. Po obliczeniu parametrów funkcji, obliczyć wartość parametru R².