Sztuczna inteligencja/SI Moduł 11 - Zadanie i metody regresji

From Studia Informatyczne


Spis treści

Metody regresji

W rozdziale tym zajmiemy się nieco dokładniej regresyjnym modelem aproksymacji, w którym aproksymowana funkcja ma postać

f: R^n \rightarrow R   \, (1)

W nieco ogólniejszym sformułowaniu, o regresji mówi się również wówczas gdy poszukuje się aproksymacji odwzorowania R^n \rightarrow R^m \,. Jeśli jednak nie zaznaczono inaczej, zakłada się postać funkcji (1).

Kwadratowa funkcja błędu

Załóżmy, że zbiór \Omega \, tworzą funkcje D\subseteq R^n \,, przy czym każda funkcja z \Omega \, jest całkowalna z kwadratem w D \,. Przyjmiemy typowe założenie, że w zbiorze \Omega \, obowiązuje metryka

|f_1, f_2| = \int_D [f_1({\mathbf x})-f_2({\mathbf x})]^2 d {\mathbf x}  \, (2)

czyli błąd kwadratowy. Założymy także, że próbkowy błąd aproksymacji jest błędem średniokwadratowym

e_S(f_1, f_2) = \frac{1}{N} \sum_{{\mathbf x} \in A} [f_1({\mathbf x})-f_2({\mathbf x})]^2  \, (3)

gdzie N \, jest liczebnością zbioru S \subset D \,, będącego zbiorem próbek.

Parametryczny model regresji

Załóżmy, że funkcja aproksymująca ma postać

h: R^n \times R^k \rightarrow R,  \, (4)

co oznacza, że przyjmuje ona n+k \, argumentów. Te dodatkowe k \, argumentów (dodatkowe w porównaniu z n \, argumentami funkcji z \Omega \,) stanowi zbiór parametrów funkcji aproksymującej, a sama funkcja h \, definiuje całą rodzinę funkcji G \subseteq \Omega \,, przy czym każda funkcja g \in G \, jest wyznaczona przez zbiór k \, wartości parametrów. W ten sposób, zadanie poszukiwania takiej funkcji g \in G \,, dla której błąd aproksymacji jest najmniejszy, da się sprowadzić do zadania poszukiwania k \,-elementowego wektora wartości parametrów. Innymi słowy, zadanie regresji parametrycznej jest zadaniem przeszukiwania przestrzeni R^k \, w celu znalezienia punktu, dla którego funkcja błędu przyjmuje minimum - punkt ten jest wektorem wartości parametrów funkcji h \,.

Funkcja reszt i gradient funkcji błędu

Powstaje pytanie - czy można podać warunki, które powinien spełniać punkt - wektor parametrów modelu - dla którego błąd jest minimalny? Odpowiedź na nie jest twierdząca.

Załóżmy, że funkcja h \, jest różniczkowalna względem swoich parametrów. Wówczas można skorzystać z warunku koniecznego istnienia minimum funkcji w punkcie, który mówi że jeśli różniczkowalna funkcja g:R^k \rightarrow R \, przyjmuje minimum w punkcie {\mathbf x}^* \,, to dla każdego i=1,...,k \, mamy \left ( \frac{\partial f}{\partial x_i} \right ) ({\mathbf x}^*) =0 \,. Innymi słowy, gradient funkcji g \, zeruje się w punkcie {\mathbf x}^* \,.

W celu znalezienia punktu, w którym funkcja przyjmuje swoje minimum, należy poszukiwać takich punktów, dla których gradient się zeruje, a jeden z nich będzie poszukiwanym punktem.

W przypadku modelu danego funkcją h \,, interesować nas będzie gradient funkcji błędu (której minimum poszukujemy) względem wartości k \, parametrów. Oznaczmy wektor parametrów przez \mathbf p \, i wektor argumentów przez \mathbf x \,, zatem funkcja aproksymująca jest oznaczana jako h({\mathbf x},{\mathbf p}) \,. Dla wygody wprowadźmy także funkcję reszt r:R^n \times R^k \rightarrow R \,, która jest dana wzorem r({\mathbf x},{\mathbf p})=h({\mathbf x},{\mathbf p}) - f({\mathbf x}) \,. Błąd średniokwadratowy dla zbioru S \subset D \, wynosi

e_S({\mathbf p})= \frac{1}{N} \sum_{{\mathbf x} \in S} [r({\mathbf x},{\mathbf p})]^2  \, (5)

Zauważmy, że błąd ten jest funkcją wartości parametrów \mathbf p \,.

Gradient próbkowego błędu średniokwadratowego względem parametrów wynosi

\frac{\partial e_S({\mathbf p})}{\partial p_i}= \frac{1}{N} \sum_{{\mathbf x} \in S} 2 r({\mathbf x},{\mathbf p}) \frac{\partial r({\mathbf x},{\mathbf p}) }{\partial p_i}  \, (6)

co, gdy uwzględnimy postać \frac{\partial r({\mathbf x},{\mathbf p}) }{\partial p_i} \,, prowadzi do

\frac{\partial e_S({\mathbf p}) }{\partial p_i}= \frac{2}{N} \sum_{{\mathbf x} \in S} r({\mathbf x},{\mathbf p}) \frac{\partial h({\mathbf x},{\mathbf p}) }{\partial p_i}  \, (7)

Reguła delta

Reguła delta jest przykładem wykorzystania warunku koniecznego istnienia minimum funkcji w punkcie, do skonstruowania prostego algorytmu poszukiwania. Zasada metody jest następująca - należy wybrać pewien początkowy wektor parametrów, obliczyć gradient funkcji błędu względem tego wektora i jeśli gradient nie jest zerowy, należy zmodyfikować wartości parametrów dodając do nich wartości wyprowadzone z gradientu. Po wykonaniu modyfikacji, powyższe czynności (począwszy od obliczania gradientu) należy powtarzać aż do osiągnięcia zerowego gradientu.

Tytułowa „delta” w nazwie reguły jest wyrażeniem

\Delta({\mathbf x},{\mathbf p}) = \frac{\partial h({\mathbf x},{\mathbf p}) }{\partial p_i}  \, (8)

i jak łatwo się przekonać, wzór na gradient błędu próbkowego można przepisać jako

\frac{\partial e_S({\mathbf p}) }{\partial p_i}= \frac{2}{N} \sum_{{\mathbf x} \in S} r({\mathbf x},{\mathbf p}) \Delta({\mathbf x},{\mathbf p})  \, (9)

Algorytm uczenia z nauczycielem, wykorzystujący regułę delta, jest podany poniżej. W algorytmie tym T \, oznacza zbiór trenujący, ||\cdot|| \, oznacza normę wektora (typowo jest to norma euklidesowa), zaś \varepsilon \, i \eta \, są parametrami liczbowymi określanym przez użytkownika, o wartościach dodatnich.

  1. Wybierz wektor parametrów {\mathbf p}^0 \,, t=0 \,.
  2. Oblicz {\mathbf d}^t= \frac{1}{N} \sum_{{\mathbf x} \in T} r({\mathbf x},{\mathbf p}^{t}) \Delta({\mathbf x},{\mathbf p}^t) \,.
  3. Jeżeli ||{\mathbf d}^t||<\varepsilon \, lub zakończ działanie, w przeciwnym przypadku przejdź do 4.
  4. Zmodyfikuj wektor parametrów {\mathbf p}^{t+1}={\mathbf p}^{t}-\eta {\mathbf d}^{t} \,
  5. Zwiększ indeks iteracji t:=t+1 \, i przejdź do 2

Uczenie polega na modyfikacjach wektora parametrów. Informacja przekazywana przez nauczyciela (który dostarcza zbioru trenującego) jest wykorzystywana do obliczenia wektora poprawek {\mathbf d}^t \,, a faza uczenia jest związana z wykonaniem kroku 4, będącego modyfikacją wektora parametrów o wektor poprawek.

Można udowodnić, że jeśli współczynnik uczenia \eta \, dąży do zera, to algorytm reguły delta jest zbieżny do punktu, którego odległość od jednego z minimów funkcji błędu dąży do zera. Potocznie mówiąc, algorytm reguły delta jest zbieżny do jednego z minimów lokalnych funkcji błędu. Jeśli minimów takich jest wiele i jeśli algorytm reguły delta jest wielokrotnie uruchamiany z różnym startowym zestawem parametrów {\mathbf p}^0 \,, to należy liczyć się z tym, że za każdym razem może być osiągane inne minimum lokalne. Przy tym warto zauważyć, że dla konkretnego startowego zestawu parametrów będzie to zawsze to samo minimum lokalne.

Model liniowy i modele „składane”

Szczególne miejsce wśród modeli regresji zajmują modele liniowe, które mają formę równania liniowego

g({\mathbf x},{\mathbf p})= \sum_{i=1}^n p_i x_i + p_0  \, (10)

gdzie {\mathbf p} \in R^{n+1} \, jest wektorem parametrów. Regresja liniowa jest omawiana zazwyczaj w czasie kursu statystyki, wspomnimy więc tylko w tym miejscu, że swoją szczególną popularność model liniowy zawdzięcza między innymi temu, że funkcja błędu ma dokładnie jedno minimum. Wystarczy więc wyznaczyć jeden zbiór wartości parametrów, dla którego gradient funkcji błędu się zeruje i tym samym mieć gwarancję znalezienia tego minimum. Przy tym metoda reguły delta nie jest najbardziej efektywna w przypadku modelu liniowego, zamiast niej stosuje się raczej metodę najmniejszych kwadratów.

Prostota modelu liniowego, stanowiąca jego niewątpliwą zaletę, nie pozwala jednak na aproksymację funkcji nieliniowych z dowolnie małym błędem. Model liniowy był jednak inspiracją do powstania całej grupy modeli nieliniowych, które nazwiemy modelami „składanymi”, gdyż każdy z nich jest liniową kombinacją funkcji nieliniowych. Dokładniej mówiąc, modele składane mają postać

g({\mathbf x}, {\mathbf P})=\sum_{i=1}^n p_{0i} \phi({\mathbf x}, {\mathbf p_i}) + p_{00}  \, (11)

gdzie {\mathbf P} \, jest macierzą (n+1)\times(n+1) \,-elementową, \mathbf p_i \, jej i \,-tym wektorem wierszowym, zaś \phi: R^n \times R^{n+1} \rightarrow R \, oznacza pewną funkcję, zwaną funkcją bazową.

W zależności od przyjętej klasy funkcji bazowych, mówi się o różnych rodzajach modeli „składanych” Poniżej wymienimy tylko niektóre z nich.

W perceptronie dwuwarstwowym funkcja \phi \, jest dana wzorem

\phi({\mathbf x}, {\mathbf p_i})=\psi\left(\sum_{j=1}^n (x_j p_{ij})+p_{i0}\right)  \, (12)

gdzie \psi: R \rightarrow R \, jest funkcją sigmoidalną - monotonicznie rosnącą, przyjmującą asymptoty skończone w minus i plus nieskończoności. Przykładami takiej funkcji mogą być funkcje \psi(x)=tanh(x) \, (asymptoty -1 \, oraz 1 \,) oraz \psi(x)= \frac{\exp(-x)}{\exp(-x)+1} \, (asymptoty 0 \, oraz 1 \,). Szerszy opis perceptronu jest podany w rozdziale 12.

W modelu RBF (ang. Radial Basis Functions - funkcje bazowe o symetrii radialnej) i w modelu falkowym (ang. wavelet) funkcja \phi \, jest dana wzorem

\phi({\mathbf x}, {\mathbf p_i})=\psi\left(\sqrt{\sum_{j=1}^n (x_j - p_{ij})^2}+p_{i0}\right)  \, (13)

Modele te różnią się wzorem funkcji \psi: R_{0+}\rightarrow R \,. W modelu RBF [1] przyjmuje się, że \psi \, jest funkcją monotonicznie malejącą do zera na dodatniej półosi, np. funkcją Gaussa [2]. Z kolei model falkowy [3] jako funkcję \psi \, przyjmuje falkę, czyli funkcję, która nie jest monotoniczna, natomiast da się ograniczyć na dodatniej półosi od góry funkcją monotocznie malejącą do zera, a od dołu funkcją monotonicznie rosnącą do zera.

Właściwość uniwersalnej aproksymacji

Powodem, dla którego interesujemy się modelami „składanymi”, jest to, że mają one właściwość uniwersalnej aproksymacji.

Definicja Właściwość uniwersalnej aproksymacji

Zbiór ciągłych funkcji G \, postaci R^n \rightarrow R \, ma właściwość uniwersalnej aproksymacji, jeśli dla każdego zwartego podzbioru K \in G \,, dla dowolnej funkcji f : K \rightarrow R \, i dla dowolnej wartości \varepsilon >0 \, istnieje funkcja g \in G \, taka, że \sup_{{\mathbf x} \in K} |f({\mathbf x}) - g({\mathbf x})| < \varepsilon \,

Innymi słowy, właściwość uniwersalnej aproksymacji zbioru G \, oznacza, że znajduje się w nim funkcja dowolnie „bliska” funkcji aproksymowanej.

Zauważmy, że nie sposób wykorzystać powyższej właściwości do konstrukcji aproksymatora. Istotnie, wprawdzie wszystkie z wymienionych powyżej modeli „składanych” mają właściwość uniwersalnej aproksymacji, lecz warunkiem jest to, żeby liczba funkcji składowych (a zatem rozmiar macierzy \mathbf P \,) była dowolnie duża. Nie jest to bardzo praktyczne rozwiązanie, zarówno ze względu na możliwość przetwarzania, jak i co ważniejsze, ze względu na niebezpieczeństwo przeuczenia rosnące wraz z rozmiarem macierzy \mathbf P \,. Nieco bardziej praktyczne wyniki niesie wymiar Wapnika-Czerwonenkisa, który, mówiąc w dużym uproszczeniu, jest wartością określającą liczbę stopni swobody modelu w zależności od liczby dostępnych danych uczących, tak aby rodzina funkcji G \, była dostatecznie bogata aby osiągnąć mały błąd aproksymacji i jednocześnie zredukować niebezpieczeństwo przeuczenia.

Przykład zastosowania - metody regresji w prognozowaniu szeregów czasowych

Szczególne miejsce z punktu widzenia częstości zastosowań metod regresji zajmują szeregi czasowe. Szeregiem czasowym jest ciąg liczb x_1, x_2,x_3,.... \, wygenerowanych jako wartości zmiennych losowych \xi_1, \xi_2, \xi_3,... \,. Zmienne losowe \xi_t \, stanowią ciąg zwany procesem stochastycznym. Nie wnikając w złożoną i ciekawą teorię procesów stochastycznych i szeregów czasowych, przyjmiemy, że między kolejnymi wartościami x_t \, istnieją pewne współzależności, a także że charakter tych współzależności nie zmienia się wraz z t \, (tzn. szereg czasowy jest stacjonarny).

Modelując stacjonarny szereg czasowy zakłada się często, że istnieje pewna funkcja f:R^n \rightarrow R \,, która pozwala na opisanie szeregu czasowego przez równanie

x_i= f(x_{i-1},...,x_{i-n})+\varepsilon_i   \, (14)

gdzie \varepsilon_i \, oznacza wartość zmiennej losowej o zerowej wartości oczekiwanej (typowo jest to zmienna losowa o rozkładzie normalnym). Tak więc przeszłe obserwacje szeregu pozwalają na określenie wartości oczekiwanej zmiennej losowej opisującej wartość przyszłą - innymi słowy, na podstawie przeszłych wartość można przewidzieć przyszłość z błędem opisanym przez rozkład \varepsilon_i \,.

Poczyniwszy takie założenie, można zadanie prognozowania wartości przyszłych szeregu czasowego przedstawić jako zadanie regresji, w którym aproksymowana jest funkcja f \, opisująca wartość oczekiwaną rozkładu przyszłych wartości. Oznacza to, że jeśli dysponujemy zapiskami historii szeregu czasowego z chwil 1,...,h \,, to jest to równoznaczne z posiadaniem h-n \, obserwacji argumentów i wartości funkcji f \,, którymi są pary:

([x_1,...,x_n], x_{n+1}), ([x_2,...,x_{n+1}], x_{n+2}),...,([x_{h-n},...,x_{h-1}], x_{h})  \,

Chcąc wykonać modelowanie szeregu czasowego, z tego zbioru obserwacji należy wydzielić zbiór trenujący i testowy. Warto zauważyć, że po przygotowaniu zbioru trenującego i testowego można proces modelowania prowadzić tak, jak w przypadku dowolnego zadania regresji, gdyż nie ma potrzeby brania pod uwagę faktu, że obserwacje zostały utworzone na podstawie szeregu czasowego.

Szereg czasowy (14) jest nazywany szeregiem autoregresyjnym (i oznaczany AR, z ang. AutoRegressive), gdyż zakłada się w jego modelu, że przyszłe wartości szeregu zależą tylko od jego przeszłości i nie mają na niego wpływu inne czynniki, w szczególności inne szeregi czasowe. Tak nie musi być. Typowym przykładem jest szereg czasowy opisujący zmiany zapotrzebowania na energię elektryczną w wielkim mieście, gdzie na typowe zachowania odbiorców energii (takich jak wieczorne zwiększenie poboru mocy przez gospodarstwa domowe) związane z trybem życia, nakładają się zależności od temperatury zewnętrznej. Pomiary temperatury zewnętrznej można również traktować jako szereg czasowy, którego jednak nie musimy modelować (chyba że chcielibyśmy w ten sposób uzyskać narzędzie prognozowania temperatury), lecz tylko uwzględniać wpływ jego przeszłych obserwacji na szereg zapotrzebowania na energię. Dostajemy wówczas model autoregresyjno-regresyjny (oznaczany ARX, z ang. AutoRegressive with auXiliary input)

x_i= f(x_{i-1},...,x_{i-na}, y_{i-1},...,y_{i-nb})+\varepsilon_i  \, (15)

w którym funkcja f \, jest odwzorowaniem f:R^{na + nb} \rightarrow R \,. Obserwacje będą miały postać par

([x_{i-na},...,x_{i-1}, y_{i-nb},...,y_{i-1}], x_{i}).  \,

Pozostałe uwagi dotyczące procesu modelowania są identyczne z poczynionymi dla szeregu autoregresyjnego. Warto również zauważyć, że w sposób naturalny można w modelu ARX uwzględniać więcej niż jeden dodatkowy szereg czasowy, w sposób analogiczny do tego, jak uwzględniono szereg y \,.

Przykład Rozważmy szereg czasowy

202, 186, 146, 138, 172, 274, 310, 434, 442, 370, 197, 259, 283, 343, 349, 383, 412, 411, 427, 378, 327, 276, 196, 166, 171, 137, 156, 185, 155, 431, 533, 815, 792, 674, 401, 429, 534, 649, 635, 671, 748, 704, 765, 671, 552, 430, 341, 262, 311, 313, 316, 333

Załóżmy, że szeregu czasowy będzie modelowany za pomocą funkcji f:R^3 \rightarrow R \,. Oznacza to, że zbiór obserwacji wartości szeregu jest w istocie zbiorem następujących par

([202,186,146],138), ([186,146,138],172), ([146,138,172],274), ([138,172,274],310), ([172,274,310],434), ([274,310,434],442), ([310,434,442],370), ([434,442,370],197), ([442,370,197],259), ([370,197,259],283), ([197,259,283],343), ([259,283,343],349), ([283,343,349],383), ([343,349,383],412), ([349,383,412],411), ([383,412,411],427), ([412,411,427],378), ([411,427,378],327), ([427,378,327],276), ([378,327,276],196), ([327,276,196],166), ([276,196,166],171), ([196,166,171],137), ([166,171,137],156),

([171,137,156],185), ([137,156,185],155), ([156,185,155],431), ([185,155,431],533), ([155,431,533],815), ([431,533,815],792), ([533,815,792],674), ([815,792,674],401), ([792,674,401],429), ([674,401,429],534), ([401,429,534],649), ([429,534,649],635), ([534,649,635],671), ([649,635,671],748), ([635,671,748],704), ([671,748,704],765), ([748,704,765],671), ([704,765,671],552), ([765,671,552],430), ([671,552,430],341), ([552,430,341],262), ([430,341,262],311), ([341,262,311],313), ([262,311,313],316), ([311,313,316],333)

Załóżmy, że pary oznaczone kolorem czerwonym są elementami zbioru trenującego, a zaznaczone kolorem zielonym - zbioru testowego.

Wykonajmy modelowanie szeregu za pomocą modelu liniowego. Model ten ma cztery parametry. W celu określenia ich wartości skorzystamy z reguły delta. Początkowe wartości przyjmijmy jako [p_0, p_1, p_2, p_3]=[0.01, 0.33, 0.33, 0.33] \,. Gradient funkcji błędu jest w modelu liniowym łatwy do obliczenia, gdyż \frac{\partial h({\mathbf x},{\mathbf p}) }{\partial p_0}=1 \, oraz \frac{\partial h({\mathbf x},{\mathbf p}) }{\partial p_i}=x_i \, dla i=1,...,n \,. W przypadku naszego przykładu gradient wynosi wynosi [-7.64, 12174.19, 6506.35, -3364.68] \,. Przyjmując współczynnik uczenia \eta=0.000001 \, otrzymamy nowe wartości parametrów [0.0100, 0.3178, 0.3235, 0.3334] \, i zaobserwujemy zmniejszenie wartości błędu średniokwadratowego. (warto tu zwrócić uwagę, że zbyt duży współczynnik uczenia, np. 0.001 \,, będzie prowadzić do zwiększenia wartości błędu). Postępując tak dalej, otrzymamy w końcu model, którego wartości parametrów prowadzą do minimalnego błędu średniokwadratowego na zbiorze trenującym. Znalezienie tego modelu oraz obliczenie wartości błędu próbkowego na zbiorze testowym pozostawiamy Czytelnikowi jako zadanie do samodzielnego wykonania.


  1. Więcej na temat modelu RBF można znaleźć w http://pl.wikipedia.org/wiki/Sieci_radialne.
  2. http://pl.wikipedia.org/wiki/Rozkład_normalny
  3. Więcej na temat modeli falkowych można przeczytać w http://pl.wikipedia.org/wiki/Falki.