Regresja liniowa wielu zmiennych (AI tutorial 2)

Poprzedni post z cyklu

Założenie, że wynik egzaminu zależy tylko od obecności na wykładach jest bardzo naiwne. Bez problemu możemy wskazać inne czynniki – ocena z ćwiczeń, IQ, nastrój prowadzącego (im lepszy tym łagodniej ocenia), ilość godzin spędzonych na nauce, etc. Uogólnijmy więc metodę regresji liniowej na wiele zmiennych. Przykłady będę prowadził na dwóch, ale podane algorytmy działają dla dowolnej liczby (w granicach mocy obliczeniowej).

Powróćmy do przykładu z egzaminem. Tym razem do dyspozycji mamy również ocenę z laboratoriów. Dostępne dane treningowe prezentują się następująco:

training

Hipoteza h(x)

Trend wydaje się następujący: im więcej nieobecności tym gorszy wynik, im lepsza ocena tym więcej punktów.
Metoda postępowania jest bardzo podobna do poprzedniej, właściwie jedyną rzeczą którą będziemy robić będzie uogólnianie. Na początek przyjrzyjmy się funkcji h_\theta(x) . Poprzednio miała postać:

h_\theta(x) = \theta_0 + \theta_1 x

Oznaczmy opuszczone wykłady przez x_1 a ocenę z laborek przez x_2. W ogólnym przypadku x_k oznacza k-tą zmienną wejściową. Funkcja będzie wyglądać następująco:

h_\theta(x) = \theta_0 + \theta_1 x_1 + \theta_2 x_2

W przypadku ogólnym (dla n zmiennych) będzie to:

h_\theta(x) = \theta_0 + \theta_1 x_1 + \theta_2 x_2 + \dots + \theta_n x_n

Spróbójmy uprościć nieco zapis (nie mam tu na myśli wstawienia znaku \sum ). Zmiennych wejściowych jest wiele, zapiszmy je jako wektor.

\boldsymbol{X} = \left [ x_1, x_2 \dots x_n \right ] n-elementowy wektor wierszowy \boldsymbol{X} \in \mathbb{R}^n

Parametrów θ jest również wiele – zróbmy z nich wektor \boldsymbol{\theta} \in \mathbb{R}^{n+1}. n + 1 ponieważ mamy \theta_0 nie powiązane z żadnym x.

\boldsymbol{\theta} = \left [ \theta_0, \theta_1 \dots \theta_n \right ]

Na razie nie pomogło nam to zbyt mocno, ilosć elementów jest różna. Wystarczy jeden prosty trick aby ją unormować. Dodajmy „fałszywą” zmienną x_0 = 1 i przypiszmy ją do pary razem z \theta_0

h_\theta(x) = \theta_0 x_0 + \theta_1 x_1 + \theta_2 x_2 + \dots + \theta_n x_n

Mnożenie przez jeden oczywiście nic nie zmieni, zostaje to co było ale w wygodniejszej formie. X wygląda teraz tak:

\boldsymbol{X} = \left [ x_0 ,x_1, x_2 \dots x_n \right ], \quad \boldsymbol{X} \in \mathbb{R}^{n+1}, x_0 = 1

Dochodzimy do konkluzji. Przypominając sobie coś algebry całą funkcje zapisujemy w postaci:

h_\theta(\boldsymbol{X}) = \boldsymbol{\theta}^T \boldsymbol{X}

Zapis matematyczny jest zwięzły, ale czy przekłada się to na kod? Tak (zakładając że korzystamy z sensownego środowiska). Używam Pythona w wersji 3.4 oraz NumPy, funkcja h prezentuje się następująco:

def h(Theta, X):
    return Theta.transpose().dot(X)

Gdzie Theta oraz X są wektorami wierszowymi. transpose() zwraca transpozycję macierzy (pamiętajmy, że wektor to szczególny przypadek macierzy gdzie jeden z wymiarów jest równy 1). dot() dokonuje mnożenia macierzowego.

Funkcja kosztu J

Nie wiem czy funkcja kosztu to poprawne tłumaczenie ale uważam je za całkiem zręczne – koszt podobnie jak funkcję J chcemy zminimalizować.

Co z funkcją J? Nic. Mamy jedynie kosmetyczną zmianę, zamiast dwóch parametrów \theta_0 oraz \theta_1 mamy teraz jeden – wektor \boldsymbol{\theta}.

J(\boldsymbol{\theta}) = \frac{1}{m} \sum\limits_{i=1}^m \left ( h_\theta(x^{(i)}) - y^{(i)} \right )^2

Również pochodne cząstkowe są praktycznie identyczne. Kontynuując nasz przykład dla dwóch zmiennych mamy trzy pochodne kierunkowe (pamiętajmy że \boldsymbol{\theta} to wektor (n+1)-elementowy gdzie n to ilość zmiennych wejściowych).

\frac{\partial}{\partial \theta_0} J(\boldsymbol{\theta}) = \frac{2}{m} \sum\limits_{i=1}^m \left (h_\theta \left (x^{(i)} \right ) - y^{(i)}\right)
\frac{\partial}{\partial \theta_1} J(\boldsymbol{\theta}) = \frac{2}{m} \sum\limits_{i=1}^m \left (h_\theta \left (x^{(i)} \right ) - y^{(i)}\right) \cdot x_1^{(i)}
\frac{\partial}{\partial \theta_2} J(\boldsymbol{\theta}) = \frac{2}{m} \sum\limits_{i=1}^m \left (h_\theta \left (x^{(i)} \right ) - y^{(i)}\right) \cdot x_2^{(i)}

Pojawiła się nowa notacja x_j^{(i)} Oznacza ona wartość zmiennej x_j w i-tym zestawie danych. Spójrzmy na przykład. Weźmy i = 8, j = 2. Na podstawie tabelki z początku postu odczytujemy że jest to 4,5 (wiersz 8, druga zmienna, dla przypomnienia wartości w kolumnie „punkty na egzaminie” to nasze y).

Pochodne po \theta_1 i \theta_2 możemy swobodnie uogólnić.
Czy to samo możemy zrobić z pochodną po \theta_0? Tak, x_0 jest przecież zawsze równe jeden. Ostatecznie wszystkie pochodne wyznaczamy ze wzoru:
\frac{\partial}{\partial \theta_j} J(\boldsymbol{\theta}) = \frac{2}{m} \sum\limits_{i=1}^m \left (h_\theta \left (x^{(i)} \right ) - y^{(i)}\right) \cdot x_j^{(i)}

Kod funkcji obliczajacej wartość pochodnej funkcji J po j-tej zmiennej dla danego \boldsymbol{\theta}

def Jpartial(j,Theta):
    sum = 0
    for i in range(m):
        sum += (h(Theta, TrainingData[i]) - Y[i]) * TrainingData[i][j]
    return 2 * sum / m

Pętlą for przechodzimy po wszystkich zestawach danych. TrainingData jest macierzą z danymi treningowymi. Dla rozważanych danych (tabela) wygląda ona tak:

TrainingData = np.array([     
	[1.,1.	,4.5],[1.,1.	,3. ],[1.,1.	,4. ],
	[1.,2.	,3.5],[1.,2.	,5. ],[1.,2.	,3. ],
	[1.,3.	,3. ],[1.,3.	,4.5],[1.,3.	,4. ],
	[1.,4.	,3. ],[1.,4.	,3.5],[1.,4.	,5. ],
	[1.,5.	,4. ],[1.,5.	,3. ],[1.,5.	,3.5],
	[1.,6.	,4.5],[1.,6.	,4.5],[1.,6.	,4.5],
	[1.,7.	,4. ],[1.,7.	,3. ],[1.,7.	,5. ],
	[1.,8.	,3. ],[1.,8.	,4. ],[1.,8.	,3.5],
	[1.,9.	,3. ],[1.,9.	,4. ],[1.,10.	,4.5],
	[1.,10.	,5. ],[1.,10.	,3. ]])
m,n = TrainingData.shape

Jest to macierz (w terminologii NumPy tablica) o m wierszach i n + 1 kolumnach. Kolumna 0 to oczywiście sztuczne x_0.
Zapis

TrainingData[i][j]

Oznacza wybranie wartości z i-tego wiersza i j-tej kolumny.

TrainingData[i]

to wybór całego i-tego wiersza.
Y to wektor z wynikami (trzecia kolumna w tabelce).

Metoda gradientu prostego dla wielu zmiennych

Już poprzednio szukanie minimum wymagało modyfikowania dwóch zmiennych. Dwie pętle for i jesteśmy w domu. Zanim to zrobimy zastanówmy się nad warunkiem stopu. Poprzednio podałem bardzo ogólne „skończ kiedy uznasz za stosowne”. Teraz zastanówmy się kiedy stosownie będzie zakończyć.

 

Analiza matematyczna mówi nam, że w optimum lokalnym wartości pochodnych są równe zero. Szansa na to, że dostaniemy *dokładnie* zero jest niewielka (a właściwie to zerowa – im bliżej jesteśmy tym wolniej idziemy, do zera zbliżamy się asymptotycznie). Pamiętajmy że obliczenia numeryczne są mniej lub bardziej przybliżone. Wybierzmy zatem pewną liczbę ε poniżej której uznamy, że funkcja się zbiegła. Jeżeli wszystkie pochodne są mniejsze od niej stwierdzamy osiągnięcie minimum i kończymy.

converge

Zbieżność do zera. X – iteracje, Y = J(θ)

notcoverage

Rozbieżność. Wartość α jest zbyt duża. X – iteracje, Y = J(θ)

 

def gradient_descent(startTheta):
    theta = startTheta
    converge = False
    i = 0
    while not converge:
        differentials = np.zeros(n)
        for j in range(n):
            differentials[j] = Jpartial(j, theta)
        for j in range(n):
            theta[j] -= alpha * differentials[j]
			
        if i % 500 == 0:
            print (np.amax(differentials))
        converge = np.amax(differentials) < epsylon
    return theta

startTheta to punkt początkowy. Dopóki nie stwierdzimy zbieżności obliczamy wartości wszystkich pochodnych. np.zeros(n) tworzy tablicę w której będziemy zapisywać ich wartości. Po ich obliczeniu przesuwamy wszystkie θ. Jeśli wszystkie są mniejsze od ε kończymy.
Dodatkowo co 500 iteracji wypisujemy informację która pozwoli nam ocenić jak długo musimy jeszcze czekać.

Wyniki

Zobaczmy jakich wartości θ wyuczył się algorytm.
\theta_0 = -0,173
\theta_1 = -0,396
\theta_2 = 1,935

Opuszczone wykłady są skorelowane z modyfikatorem -0,4. Im więcej opuszczonych tym mniej punktów. Ocena z laborek ma mnożnik ~2. Widać że komputer doszedł do „rozsądnych” wniosków (na podstawie zadanych wartości treningowych). Zauważyć można również fakt, że ocena jest ~4x bardziej istotna przy przewidywaniu wyniku. Zapytajmy teraz o kilku studentów i sprawdźmy co na temat ich egzaminu powie algorytm.

results

Pamiętajmy, że jakość otrzymywanych wyników zależy w kluczowym stopniu od podanego zestawu treningowego!

Całość kodu

Reklamy

2 comments

Skomentuj

Wprowadź swoje dane lub kliknij jedną z tych ikon, aby się zalogować:

Logo WordPress.com

Komentujesz korzystając z konta WordPress.com. Wyloguj / Zmień )

Zdjęcie z Twittera

Komentujesz korzystając z konta Twitter. Wyloguj / Zmień )

Zdjęcie na Facebooku

Komentujesz korzystając z konta Facebook. Wyloguj / Zmień )

Zdjęcie na Google+

Komentujesz korzystając z konta Google+. Wyloguj / Zmień )

Connecting to %s