Programowanie

Regresja logistyczna (logistic regression) – AI tutorial 3

Pierwszy post z cyklu
Drugi post z cyklu

Zostawmy na jakiś czas egzaminy, przejdźmy do innych problemów.
Idąc do lekarza mówimy co nas boli, robimy badania (jeśli NFZ przypadkiem jeszcze nie zbankrutował) a wyniku pracy inteligentnej osoby dostajemy odpowiedź czy jesteśmy chorzy. Załóżmy, że chcemy zastąpić lekarzy maszynami. Doktor na podstawie zdobytej wiedzy oraz doświadczenia (czyli de facto również wiedzy) podejmuje decyzję, czy jesteśmy zdrowi czy chorzy. Jak to się ma do poprzedniego problemu? Trochę podobnie, trochę inaczej. Na wejściu jest to samo – n zmiennych wejściowych (po ang. zwane features). Wyjście jest inne – zamiast dostać wartość z ciągłego zbioru (punkty) dostajemy binarną odpowiedź chory/zdrowy. Mamy tutaj do czynienia z problemem klasyfikacji – efekt pracy może należeć tylko do pewnego skończonego zbioru. Wymaga to nieco innego podejścia. Inne przykłady takich problemów:

X oznacza wektor zmiennych wejściowych

  • Email klasyfikowany jako spam X → {0, 1}
  • Email klasyfikowany do pewnej kategorii (praca, prywatne, newslettery) X → {0, 1, 2}
  • Kategoria aut (sedan, pickup etc) na podstawie zdjęcia X → {0, 1, 2, 3, 4…}
  • Klasyfikacja płci X → {0, 1}

Testowany obiekt musi być w dokładnie jednej kategorii.

Na razie ograniczmy się do odpowiedzi tak/nie: y \in \{0,1\}. Postawmy taki problem: mamy chorego na określoną chorobę i chcemy określić czy może zostać wyleczony czy też nie. Do dyspozycji mamy takie dane jak jego wiek, waga, wyniki badań, płeć etc. Na początku załóżmy że rozważamy tylko wiek. Przykładowe dane mogą wyglądać następująco:

pic1

Kółko oznacza 0, gwiazdka 1 – będę dalej używać tej notacji.

Widzimy, że tendencja jest bardzo prosta – im pacjent starszy tym gorsze rokowania. Gdybyśmy chcieli wyznaczyć przedziały dla y=0 (wyzdrowienie) oraz y=1 (śmierć) mogło by wyglądać to tak:

pic2

Widać tutaj pewien problem – w okolicach miejsca podziału pewnie obiekty zostaną źle sklasyfikowane, nie da się tego całkowicie uniknąć. Co możemy z tym zrobić? Zamiast mówić „Hej, umrzesz” lepiej będzie powiedzieć „Umrzesz z prawdopodobieństwem 87%”. Oczywistą zaletą jest, że w niepewnej okolicy będziemy dostawać wyniki oscylujące wokół 50% co wskaże fakt niepewności klasyfikacji. Rzućmy okiem na przykład z dwiema zmiennymi wejściowymi:

db1

Jedna zmienna jest na osi X, druga na Y. Zaznaczona linia wskazuje miejsce w którym tniemy płaszczyznę (czy to zawsze i wszędzie muszą być linie proste? Nie, pokażemy wkrótce inne warianty).

Hipoteza h(x)

Ostatnio szukaliśmy takiej funkcji h która dla pewnego zestawu danych zwróci nam przewidywany wynik egzaminu. Teraz szukamy takiej która zwróci nam prawdopodobieństwo tego, że y = 1. W związku z tym funkcja musi być ograniczona 0 \le h_\theta(x) \le 1 . Konstrukcja takiej funkcji może być bardzo problematyczna.
Prościej było by uznać że jeśli h(x) > 0 – wynik pozytywny (y = 1). W takim rozumieniu tracimy jednak informację o prawdopodobieństwie chociaż intuicyjnie czujemy że im większa wartość tym większa pewność wyniku pozytywnego, a im mniejsza negatywnego. Potrzebujemy więc pewnej funkcji (nazwijmy ją S) która zamieni wygodną dla nas reprezentację na tą bardziej użyteczną.

S: (-\infty; +\infty) \mapsto (0;1)

Taka funkcja istnieje, nazywa się ona sigmoid function (tłumaczenia się nie podejmę) i jest zdefiniowana następująco:

S(t) = \frac{1}{1+e^{-t}}

Spójrzmy na wykres

 

Nada się znakomicie – dla t = 0 zwraca dokładnie 0,5. im dalej w kierunku nieskończoności tym bliżej 0 lub 1.

Uporządkujmy nieco oznaczenia. Funkcja hipotezy h_\theta(x) będzie teraz korzystać z funkcji S aby „przetłumaczyć” wynik z formy wygodnej dla nas na tą użyteczną. Do kreślenia prostych używaliśmy ostatnio h_\theta(\boldsymbol{X}) = \boldsymbol{\theta}^T \boldsymbol{X} . Teraz też będzie odpowiednia. Oznaczmy ją tym razem jako g(X):
g(\boldsymbol {X}) = \boldsymbol{\theta}^T \boldsymbol{X}
pamiętając że mamy do dyspozycji funkcję S:
S(t) = \frac{1}{1+e^{-t}}

Otrzymujemy ostatecznie:

h_\theta(\boldsymbol{X}) = S(g(\boldsymbol{X}))
lub równoważnie:
h_\theta(\boldsymbol{X}) = \frac{1}{1+e^{-\boldsymbol{\theta}^T \boldsymbol{X}}}

Kiedy będziemy chcieli zmienić zachowanie funkcji zmodyfikujemy tylko g(X)

Właśnie, jak mamy interpretować wartość zwracane przez g? Przy regresji liniowej była to wartość oczekiwana dla danego wejścia (co zresztą było bardzo dobrze widać na wykresie). Do tej pory ustaliliśmy: funkcja h zwraca prawdopodobieństwo wyniku pozytywnego (y = 1). Jeżeli jest ono większe niż 50% klasyfikujemy wynik jako pozytywny, w przeciwnym razie jako negatywny. Progiem „przełączania” jest wartość 0.5 funkcji h czyli wartość 0 funkcji g. Wynika z tego, że równanie g(x) = 0 pokazuje linię (powierzchnię, hiperpłaszczyznę etc) po której jednej stronie wartości są klasyfikowane jako pozytywne, a po drugiej – negatywne. Spójrzmy na przykłady

db2

 

Dotychczas używana przez nas metoda (konkretnie funkcja g) pozwalała tylko na oddzielanie elementów linią prostą. Nic jednak nie stoi na przeszkodzie aby była nieco bardziej skomplikowana, np:

g(x) = \theta_0 x_0 +\theta_1 x_1+\theta_2 x_2+ \theta_3 x_1^2 +\theta_4 x_2^2 +\theta_5 x_1 x_2 + \dots

Spowoduje to dopasowanie pewną krzywą, im dłuższe wyrażenie tym potencjalnie bardziej skomplikowaną. Rozwiniemy temat przy okazji sieci neuronowych (przykład na coś innego niż prosta pojawi się jeszcze w tym poście).

Funkcja kosztu J

Funkcję hipotezy mamy za sobą, co z funkcją kosztu J? Poprzednio miała formę:

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

Wyizolujmy z niej część odpowiedzialną za obliczenie kosztu pojedynczego zestawu danych. Będzie to:
h_\theta(x^{(i)}) - y^{(i)}

Zastanówmy się jak powinna wyglądać obecnie. Z całą pewnością chcemy, aby poprawny wynik był obarczony zerowym kosztem. Z drugiej strony jeżeli nasza hipoteza stwierdzi z dużą pewnością nieprawidłowy wynik, penalizacja powinna być wysoka. Pamiętajmy że poruszamy się w przedziale (0;1).

Na początek załóżmy, że y = 0. Jak ma kształtować się kara dla różnych wyników h(x)? Na przykład tak:

j0

Dla prawidłowego rezultatu kara jest równa zero. Im bardziej algorytm był pewny złej odpowiedzi tym większą karę ponosi (y dąży do nieskończoności wraz z x dążącym do 1).

Dla y = 1 sytuacja jest odwrotna:

j1

Wyrażenie opisujące pojedynczy zestaw danych będzie więc wyglądać tak:

\text{single penalty}(\boldsymbol{X}) =  \begin{cases}  -\log(h_\theta(\boldsymbol{X})) & \text{if } y = 1 \\  -\log(1-h_\theta(\boldsymbol{X})) & \text{if } y = 0  \end{cases}

Możemy pozbyć się niepożądanego warunku zapisując całość jako: (ze względu na y \in \{0,1\})

\text{single penalty}(\boldsymbol{X}) =-y \log(h_\theta(\boldsymbol{X})) - (1-y) \log(1-h_\theta(\boldsymbol{X}))

Podstawiając to do poprzedniego wyrażenia na wartość funkcji J(\boldsymbol{\theta}) otrzymujemy ostatecznie:

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

Optymalizacja

Tradycyjnie, kiedy mamy już opisną funkcję J problem jest sprowadzony do jej minimalizacji. Wyznaczmy pochodne kierunkowe względem wektora \boldsymbol{\theta}.

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

???

Nie, to nie jest błąd. Czy dostaliśmy to samo (z dokładnością do mnożenia przez 2)? Prawie. Pamiętajmy że obecne h_\theta(x) wygląda inaczej niż w regresji liniowej.

Całość algorytmu gradientu prostego jest identyczna, więc nie będę go opisywał ponownie.

Implementacja

Zaimplementujmy wersję w której funkcja g(x) będzie prezentować się następująco:

g(x) = \theta_0 x_0 + \theta_1 x_1^2 + \theta_2 x_2^2

Jest to nic innego jak elipsa o środku w O=(0,0).

Funkcje S oraz h w kodzie Pythona:

def sigmoid(x):
    return 1/(1 + np.e**(-x))

def h(Theta, X):
    z = Theta[0] * X[0] + Theta[1] * X[1] **2 + Theta[2] * X[2]**2
    return sigmoid(z)

Obliczenia pochodnych cząstkowych są identyczne jak poprzednio podobnie jak i sam gradient. Spójrzmy na przykładowy wynik działania. Mając dostępne następujące dane:

raw

Algorytm wytyczył krzywą dla której prawdopodobieństwo obu wyników jest takie samo (P(y=1)=0,5 w sposób następujący:

line

Dodatkowo krzywe dla 75% prawdopodobieństwa y=1 oraz y=0.

imgNastępnym razem zostawimy uczenie się na pewien czas i zastanowimy się nad podejmowaniem rozsądnych decyzji.

Całość kodu

Reklamy

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

It’s been a long time (AI tutorial 1)

Sporo czasu upłynęło od ostatniego razu. Po raz kolejny wznawiam bloga, tym razem mając konkretne cele. Blog ma mi towarzyszyć w moich zmaganiach na drodze do zostania wymiataczem w AI. Temat bardzo ciekawy, bardzo trudny, bardzo nęcący ale zawsze odstawiany przeze mnie na bok ze względu na różne wymówki. Wystarczy, w końcu coś pęka i człowiek zaczyna robić to co powinien zacząć już dawno. Każda kolejna rzecz opisana na blogu w celu raz usystematyzowania wiedzy, dwa posiadania notatek.

Na początek odświeżenie tego co już wiadomo.

Wyobraźmy sobie taką sytuację. Jest egzamin, jedni zdają inni nie zdają, wszystko w normie. Każdy jednak woli być w tym zbiorze studentów który zdał. A najlepiej było by wiedzieć czy się zdało jeszcze przed. Można wykorzystać dane historyczne aby spróbować przewidzieć przyszłość. Czy można to nazwać AI? IMO tak. Naturalne jest że wiele decyzji podejmujemy na postawie poprzednich doświadczeń. Niech komputer na podstawie tego co zdarzyło się wcześniej przewidzi przyszłość. Możemy to nazwać uczeniem maszynowym. Pewne dane są dostępne i zrozumiane. Przekazujemy je do programu, który na ich podstawie uczy się pewnych wzorców. Kiedy skończy podajemy nowe dane a on przekazuje przewidywany wynik. Brzmi nieźle.

Drawing1

Zacznijmy od czegoś prostego.

Mamy dostępne dane na temat zdobytych punktów z egzaminu oraz opuszczonych wykładów. Na tej podstawie spróbujmy przewidzieć jaki wynik uzyska ktoś kto był dokładnie na połowie wykładów.
Zobaczmy jak to wygląda na wykresie.

dd

Jak widać dane układają się mniej więcej wzdłuż pewnej linii. Spróbujmy ją naszkicować na oko.

łan

Widać że dla 50% obecności możemy oczekiwać wyniku w okolicach 5 punktów. No dobrze, ale jak komputer ma do tego dojść?

Regresja liniowa jednej zmiennej.

Podsumujmy co mamy do tej pory: zbiór par wartości (obecność, punkty). Szukamy pewnej funkcji liniowej która pozwoli przewidywać wynik. Postać ogólną, którą wszyscy znają wygląda tak:
y = ax + b
Dalej będę używał formy:
h_\theta(x) = \theta_0 + \theta_1 x

Potrzebujemy wyznaczyć dwa współczynniki \theta_0 i \theta_1 . Tym zajmie się tytułowy algorytm.

Aby stwierdzić czy rozwiązanie jest dobre zdefiniujmy sobie coś co pozwoli ocenić jego jakość. Taka funkcja musi sprawdzać jak bardzo rzeczywiste dane różnią się od przewidywanych. Najprościej będzie użyć sumy kwadratów różnic, innymi słowy metody najmniejszych kwadratów.

Dla pojedynczej pary wyrażenie będzie wyglądać tak:
(h_\theta(x) - y)^2
Podnosimy do kwadratu różnicę między tymi wartościami. Gratisowo pozbywamy się ew. wartości ujemnych.
Jak będzie to wyglądać dla całego zbioru danych?

\sum\limits_{i=1}^m (h_\theta(x^{(i)}) - y^{(i)})^2
Gdzie x^{(i)} i y^{(i)} oznaczają po prostu i-tego x i y, a m liczność zbioru z danymi.

Ostatecznie zapiszmy to jako średnią powyższych odchyleń:
J(\theta_0, \theta_1) = \frac{1}{m} \sum\limits_{i=1}^m \left ( h_\theta(x^{(i)}) - y^{(i)} \right )^2

Warto zauważyć, że h_\theta(x^{(i)}) to to samo co \theta_0 + \theta_1 x

Będziemy starali się tak manipulować parametrami \theta_0 i \theta_1 aby zminimalizować funkcję J i w efekcie uzyskać lepsze dopasowanie.

Minimalizacja

Pamiętajmy, że rozwiązanie problemu sprowadza się teraz do znalezienia minimum funkcji J.

Funkcja J jest dwuargumentowa. Do wyznaczenia minimum (niekoniecznie globalnego!!!) posłużymy się algorytmem gradientu prostego. Nazwa sugeruje że jest to prosta metoda i tak jest w istocie.

Pomysł jest następujący:
1. Wybierz punkt początkowy (wartości dla \theta_0 i \theta_1 ).
2. Rozejrzyj się dookoła gdzie możesz zejść w dół (i zejdź).
3. Powtarzaj tak długo aż znajdziesz optimum/skończy się czas/etc.

Krok 1 jest oczywisty. Zastanówmy się nad drugim. Aby zdecydować w którym kierunku iść, będziemy potrzebowali pochodnych cząstkowych (jak zmienia się funkcja w obecnym położeniu). Wyglądają one następująco:

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

W jaki sposób pochodne mają nam pomóc? Aby zobaczyć jak to działa weźmy funkcję kwadratową.
y=x^2-2x+3
Jej pochodna
y'=2x-2
Weźmy punkt startowy x = 10 i poszukajmy minimum. Będzie nam potrzebny jeszcze współczynnik \alpha . O nim za chwilę, przyjmijmy że ma wartość \alpha=0.1

Pochodna w tym miejscu jest równa y'(10) = 18
Zmodyfikujmy x następująco:
x_{i+1} = x_i - \alpha \cdot 18 1iter
Przesunęliśmy się we właściwym kierunku!
Wykonajmy więcej iteracji.

10iter

Jak widać osiągnęliśmy optimum. Czym jest \alpha? Jest to współczynnik mówiący jak ostro mamy iść w kierunku rozwiązania. Im mniejszy tym więcej iteracji trzeba wykonać, im większy tym szybciej ale musimy uważać na zbyt dużą wartość. Ustawmy \alpha=1.2, 20 iteracji.

overshoot

Coś nie zagrało. Nad przyczyną polecam zastanowić się samemu.

OK, sprawdziliśmy że metoda działa, przełóżmy ją na dwie zmienne.
Rozumowanie jest identyczne, jedyne o czym musimy pamiętać to jednoczesna zmiana obu parametrów (liczymy krok dla \theta_0 oraz \theta_1 a dopiero potem zmieniamy wartości tych zmiennych).

Co zobaczymy po implementacji?

lech

Seems legit

 

Zobaczmy jak będzie wyglądał kod.

X = [10,9,8,7,6,5,4,3,2,1,0]
Y = [1,2,4,3,5,4,6,7,4,8,9]

m = len(X)
alpha = 0.01

Wprowadzamy dane, ustawiamy m oraz α.

def h(x, theta):
    return theta[0] + theta[1] * x


def J_partial0(theta):
    sum = 0
    for i in range(m):
        sum += h(X[i], theta) - Y[i]

return 2 * sum / m;

def J_partial1(theta):
    sum = 0
    for i in range(m):
        sum += (h(X[i], theta) - Y[i]) * X[i]

    return 2 * sum / m;

Definiujemy fukcję h oraz pochodne funckcji J.

def gradient_descent(startTheta):
    theta = startTheta
    for i in range(5000):
        d0 = J_partial0(theta)
        d1 = J_partial1(theta)
        theta[0] -= alpha * d0
        theta[1] -= alpha * d1

    return theta

Pamiętamy aby obliczyć d0 oraz d1 przed ich dodaniem do aktualnego wyniku.

Podsumowanie

To wszystko, metoda jest bardzo prosta ale to dopiero pierwszy krok. Następnym razem: co jeśli mamy dodatkowe dane i chcemy je uwzględnić w uczeniu.

Cały kod dostępny tutaj: https://gist.github.com/Xevaquor/e128c3c31c48dafd5766

DISCLAIMER: Używana notacja jest luźna, niekoniecznie musi być w pełni formalnie poprawna.

ThunderWarn

Zainspirowany wpisem Gethioxa pomyślałem o napisaniu monitora nastroju Zeusa pod Windows. Założenia są podobne – pobieranie mapki z burzami oraz jej wyświetlenie na pulpicie. I tutaj jest pierwszy poważny problem – Gethiox rysuje mapę na wyświetlanej tapecie. Chciałem osiągnąć podobny efekt ale bez bazgrania po tapecie – póki co nie mażę po niej ale podobnego efektu nie uzyskałem –  wszystko jest wyświetlane w okienku.

Integracją z pulpitem bym tego nie nazwał…

Myślałem o użyciu ActiveDesktop ale to jest już od dawna niewspierane. Gadżety pulpitu też nie wyglądają na technologię której użycie nie jest równoważne nekromancji. Póki co jest to okienko WPF. Właśnie WPF – nie wiem co mam myśleć o tym wynalazku – niby już dobry czas temu czytałem że jest to zalecana technologia tworzenia UI, w przeciwieństwie do niezalecanego WindowsForms – ale nie mogę oprzeć się wrażeniu że sam MS nie jest za bardzo tym zainteresowany i nie bardzo wiadomo co tym wynalazkiem zrobić (przenieść do Silverlight który miał być rewolucją a potem go zabić, a XAML przepchać dalej do Metro (czy tam Modern UI jak to się teraz nazywa)). Co do samego WPF – większość czasu który spędziłem do tej pory polegało na rozgryzaniu o co w tym bałaganie chodzi. W moim odczuciu jest to dość dziwna technologia, gdzie zrobienie prostych rzeczy jest zdecydowanie zbyt skomplikowane, całość jest przekombinowana a wczytanie obrazka do pamięci (do pamięci, nie trzymając blokady na pliku) jest co najmniej nietrywialne. O banalnym przeniesieniu danych z System.Drawing.Bitmap z GDI+ do tego cuda nie wspominając.

W porównaniu do oryginalnego skryptu  kod jest wykładniczo bardziej skomplikowany. Jest to spowodowane po pierwsze technologią a pod drugie docelową większą funkcjonalnością. W założeniu będą dostępne dodatkowe źródła danych poprzez wtyczki, ostrzeganie przez zbliżająca się burzą i automatyczną ustalanie lokalizacji. Przydało by się trochę piorunów do testów. W chwili obecnej odpowiada w zasadzie inspiracji (mając ~10 razy większą objętość, pomijając kod wygenerowany przez VS). Jest to fajna okazja żeby pobawić się odważniej gałęziami Gita oraz innymi rzeczami. Bardzo chciałem wykorzystać async/await ale na moim Win8 RP, Visual w wersji RTM nie ma zamiaru działać – zamiast tego wykorzystałem więc TPL do obsługi asynchroniczności (nie bez problemów). Bajerami z nowego dotnetu pobawię się kiedy w końcu nowy system będzie miał premierę.

Caly kod dostępny jest na GitHub – https://github.com/Xevaquor/ThunderWarn

Developer’s Adventure – 0x01

Wakacje się skończyły, okres wypoczynku też – czas powrócić do pracy. Praca na WSoC ma w tej chwili bezwarunkowo najwyższy priorytet – został raptem miesiąc na skończenie gry. Kilka dni temu nie było jeszcze właściwie nic – poza asynchronicznym ładowaniem zasobów (jestem z tego cholernie dumny!) i przygotowaniem kodu pod przerzucanie stanów gry. Przez ostatnie trzy dni udało się zintegrować silnik fizyczny Farseer który jest portem Box2D dla XNA, oraz zaimplementować pomniejsze funkcjonalności. Ku mojemu zaskoczeniu nie miałem problemu ze zrozumieniem swojego kodu po pięciotygodniowej przerwie – jednak potrafię uczyć się na błędach i mój kod jest coraz łatwiejszy do czytania.

Przy tej okazji miałem okazję zapoznać się z narzędziem FxCop – służy ono do statycznej analizy kodu pod względem różnego rodzaju błędów, albo potencjalnych błędów. Czasem mam wrażenie że zbytnio się czepia, ale okazało się użyteczne wskazując mi niezauważone miejsca w których mogły wystąpić problemy. Innym narzędziem jest Visual Studio w wersji Ultimate. Oczywiście 90-dniowy trial. Na razie miałem okazję bliżej zapoznać się tylko z profilerem oraz InteliTrace. Jest to tylko niewielki wycinek, a mimo to jestem pod wrażeniem potęgi Visuala. Skoro o narzędziach mowa Developer’s Adventure jest pierwszym projektem w którym ‚na poważnie’ używam Gita. Kilkukrotnie okazał się bardzo pomocny w sytuacjach typu „przecież jeszcze przed chwila to działało”. Nie wiem jak mogłem wcześniej bez tego pisać kod.

Warsztat Summer of Code

Ponad rok temu miałem definitywnie zakończyć swoją przygodę z gamedevem. Przez rok nawet nie tykałem się pisania żadnej gry – nawet gdybym chciał nie miałem na to czasu. To nie jest branża dla mnie – nadal jestem tego całkowicie pewny.

Pomimo to postanowiłem spróbować swoich sił w tegorocznym Warsztat Summer of Code. Sam nie wiem dlaczego – ale warsztat ma jakąś magię która mnie tam przyciąga i nie pozwala odejść. Żeby jednak nie odchodzić zbyt od podstawowych narzędzi pracy gra powstanie w C# i XNA (sorry Linux). Nie liczę na żaden spektakularny sukces, właściwie sam nie wiem na co liczę biorąc udział. Może będzie to odskocznia od kilku nudnych projektów które muszę zrealizować (albo raczej pretekst żeby ich nie ruszać).

Teraz słowo o grze (a raczej wklejony plik żeby nie pisać tego samego n-ty raz)

Gra pisana na Warsztat Summer of Code 2012 – zapewnie nic z tego nie wyniknie tak jak co roku ale tak jak co roku spróbuję po raz n-ty.

//Wszelkie podobieństwo do Warsztat Game niezamierzone i całkowicie przypadkowe.

Hack’n’slash w którym wcielamy się w młodego programistę ktory musi stawić czoła wielu
niebezpieczeństwom aby wygrać. Gracz niszczy potwory strzelając w nie bliżej nieokreślonymi pociskami
(no bo czym ma strzelać? pudełkami po pizzy i kubkami po kawie?). Gra podzielona na niezależne
poziomy. W nagrodę za kolejne sukcesy gracz dostaje bonusy dające mu większą moc niesienia pożogi,
ale i przeciwnicy są coraz trudniejsi (ciekawe kto to zbalansuje?). Widok 2D z góry.

Gracz:
Do wyboru dwie/trzy klasy postaci: programista C++, Java i może Assembly (można by i więcej
ale muszą się jakoś sensownie różnić). Przykładowe różnice:
-> Java – wolny i słaby w ataku ale odporny na ciosy
-> Asm – szybki i z wielką mocą ale podatny na ciosy
-> C++ – coś po środku (żeby się tylko nie okazało że jedyna grywalna postać)

Gracz wali jakimiś pociskami (dla Javy te kubki z kawą to by nawet pasowały),
każda klasa ma inny kolorek i inną moc. Po zdobyciu iluś punktów
(albo lepiej EXP – lepiej to wygląda) gracz zdobywa dodatkowe moce.

Moce:
//DRZEWKO Z MOCAMI!!!
(przykładowe wszystko)
Wątki – można walić więcej niż jeden pocisk na raz
Sieci – rażenie kilku wrogów na raz
Asercje(debugger) – leczenie/regeneracja HP
Unit testy – ochrona przed słabymi przeciwnikami (od buga w dół)

Przeciwnicy:
N00by – małe, niegroźne. Jednynie co to upierdliwe.
Trolle – Trochę groźniejsza wersja
Anoni – mocniejsza wersja powyższych.
Bugi – silny przeciwnik (it’s not bug – it’s feature <- fajnie by to jakoś wykorzystać)
Memleaki – potężniejsza wersja bugów
Project Manager – największe zuo, bardzo groźni ale nieliczni
Przeciwnicy pojawiają się na mapie w losowych miejscach na starcie + spawnują się
aby zniechęcić gracza do bierności.

Poziomy:
Każdy level ma nazwę w postaci błędu np:
stack overflow
unresolved external symbol
access violation
Generowane proceduralnie – na pewno nie będzie mi się chciało pisać edytora
Im późniejszy level tym trudniejszy (proporcje i ilość wrogów)

Inne uwagi:
Wplecenie komiksów xkcd między poziomy.

I tak się na 90% nie uda ale cóż.
Projekt jest dostępny na GitHub: https://github.com/Xevaquor/adventure

Temat na Warsztacie: http://forum.warsztat.gd/index.php?topic=25555

Pewnie po raz kolejny popełniam ten sam błąd, ale cóż 🙂

Quo vadis?

Ostatni miesiąc przeznaczyłem głównie na rozmyślania nad swoimi celami. Tutaj podzielę się tym co wymyśliłem w sferze związanej z programowaniem. W trakcie mojej krótkiej aczkolwiek burzliwej kariery zapoznałem (mniej lub bardziej) pobieżnie się z wieloma językami i zastosowaniami programowania. Były to głównie:

  • HTML + CSS + JavaScript
  • PHP/Ruby on Rails
  • Python
  • C# i .NET oraz XNA
  • Java
  • C++ i SDL/Allegro/SFML/DirectX/OpenGL/Irrlicht/OGRE
  • PowersShell i Bash
  • Assembler
  • Delphi/Object Pascal

Jak widać rozrzut technologii dosyć spory ale sądzę, że dało mi to wystarczająco dużo doświadczeń aby móc zdecydować na czym się skupić (jak coś jest do wszystkiego to jest do niczego). Epizod związany z gamedev był chyba najdłuższy i najbardziej zaawansowany, ale ostatecznie wiem że nie chcę się tym zajmować. Technologie webowe również mnie nie pociągają.

Wybrałem sobie trzy które intensywnie będę zgłębiać dalej:

  1.  Assembler – tak wiem, że tego się już nie używa, że trudny, że błędogenny, że długo się w nim pisze a rezultaty nie są spektakularne. Ale assembler ma dla mnie jedną niezwykle istotną cechę: jest na bardzo niskim poziomie abstrakcji co pozwala mi zrozumieć jak wszystko dokładnie działa. Pisząc w C++ wiedziałem, że „argumenty wrzuca się na stos od końca”. Pisząc w assemblerze wiem, że „argumenty umieszczane są na stosie instrukcją push, która… […] …rejestr esp wskazuje na ostatni element stosu…”. Po co mi to wiedzieć? Prawdę powiedziawszy nie wiem, ale czuję się dużo bardziej komfortowo rozumiejąc co dzieje się za kulisami.
  2. C++ – w tym języku programuję w sumie najdłużej i najintensywniej. Podoba mi się w nim jego względna wysokopoziomowość programowania obiektowego jak względna niskopoziomowość (wstawki asm, wskaźniki ❤ ). Ponadto w tym języku zwykle widać co się tak naprawdę dzieje, można swobodnie kontrolować czy przekazać parametr przez wartość czy przez referencje oraz stosować magiczne szablony i makra 😉
  3.  Java – rozwiązanie wysokopoziomowe. W tej kategorii wygrała dzięki swojej przenośności. Javę mają prawie wszyscy, ile znacie osób które na Windows mają Pythona? .NET to z kolei w mojej opinii takie trochę niewiadomo-co. Ani to przenośne, ani przywiązane do Windows (do użycia natywnych funkcji trzeba użyć platform invocation). Kończąc wątek .NET szkoda mi C# w którym pisze mi się doprawdy sympatycznie (lambdy, LINQ). W zasadzie dopiero zaczynam poznawać Javę lecz sugerując się powszechną opinią odnośnie jej podobieństwa do C# jestem nastawiony optymistycznie.

Z niecierpliwością czekam na nowy PC (to będzie najdłuższy miesiąc w moim życiu 😉 ). Wtedy zacznę pisać super-pure-64-bit-mega-ultra OS który na wzór i podobieństwo Mac OS będzie kompatybilny tylko sam ze sobą 😉 .

Koniec przerwy

Jak widać ostatni post jest sprzed dwóch miesięcy. Planowałem wtedy zmienić jezyk na angileski z czego jak widać zrezygnowałem. Powodem jest fakt, że przez cały czerwec nie zdołałem napisać nic w tym jezyku z czego byłbym zadowolony. W lipcu nie napisałem nic gdyż pierwszy miesiąc wakacji przeznaczyłem na ozważania metafizyczne i beztroski odpoczynek. W tym także od programowania – przez połowę wakacji nie napisałem nawet 200 linijek.
Sierpień będzie bardziej pracowity. Głownym celem jest nauczyć się wreszcie pisać bezwzrokowo (wstyd się przyznać ale nie potafię!). Dodatkowo mam zamiar ogarnąć w pewnym stopniu PowerShell (o którym zapewne wkrótce napiszę).
Teraz zamykając wątek Maasym the Mage. Nie chcę nawet o tym myśleć ale zamknę ten rozdział ostatecznie. Ostatnia wersja jest do pobrania tutaj:
http://dl.dropbox.com/u/19656258/Maasym%20the%20Mage%20v.0.88%20FINAL.rar
Źródła:
http://dl.dropbox.com/u/19656258/Maasym%20the%20Mage%20v.0.88%20FINAL%20-%20SOURCE.7z
Gra wcale nie wyglada i nie działa tak jak miało to wyglądać, ale faktem jest że organizatorzy konkursu też się nie popisali gdyż tuż przed końcem amin doszedł do wniosku „że już mu się nie chce”, forum było częściowo niedziałające, i tak naprawde czart jeden wiedział czy konkurs dojdzie do końca. Gdyby ktoś chciał „pociągnąć” projekt dalej to być może będę w stanie udzielić informacji odnośnie istniejącego kodu (na tyle na ile pamiętam).
Tyle na dzisiaj.

Maasym the Mage

Miesięczna przerwa na blogu nie jest powodem do dumy (wpisy miały być co najmniej raz na tydzień) więc muszę nadrobić zaległości 😉

Tak się złożyło, że postanowiłem wziąć udział w konkursie DevLogs. Była to jedna z moich najlepszych decyzji w ostatnim czasie. Dlaczego – o tym później. Mój projekt nazywa się Maasym the Mage i w skrócie przedstawia się jako ‚platformówka 2d w której chodzimy magiem i odgrywając melodyjki rzucamy czary zabijając przeciwników, gromadząc expa, ucząc się nowych czarów a przy okazji mamy okazję uratować świat przed zagładą. Nie ma sensu powtarzać dwa razy tego samego więc w tym miejscu podlinkuję devlog. Z bardzo dużą dozą prawdopodobieństwa mogę powiedzieć, że to będzie moja ostatnia próba napisania gry. Kończenie projektów zawsze sprawiało mi ogromne problemy, ale ten będzie ciągnięty aż do ukończenia. Dlaczego? Bo chcę ostatecznie przelać czarę goryczy.

Niestety ale, programowanie gier z całą pewnością nie jest tym co chcę robić w życiu. Szkoda tylko, że ta refleksja przyszła tak późno. No cóż, mimo wszystko lepiej teraz niż za rok. Gry były bodźcem, który pchnął mnie do programowania. Jednak nie są czymś co może mnie utrzymać przy monitorze cały dzień. Fakt, że nie zauważyłem tego wcześniej był zapewne spowodowany tym, że nie kończyłem całkowicie swoich gier. Pojawiały się różne myśli, ale brakowało wystarczająco silnego bodźca. Maasym narodzi się w bólach choćby nie wiem co. Ale dzięki niemu próg cierpienia(*) zostanie przekroczony, a ja wyzbędę się wszelkich wątpliwości.

* – to za duże słowo, ale aktualnie nie przychodzi mi na myśł żadne lepsze.

Data Breakpoint w Visual C++

Co to jest breakpoint (pułapka) wie każdy, nawet początkujący programista. Przerwanie wykonywania programu w momencie osiągnięcia pewnego miejsca jest niezwykle przydatne w trakcie debugowania – zwłaszcza kiedy program wywala się z niewiadomych powodów. Najczęściej wykorzystywane w Visual Studio są zwyczajne breakpointy wstawiane poprzez kliknięcie na lewym marginesie kodu. W dzisiejszym wpisie zajmę się drugim rodzajem dostępnym w wersji Express (2008). Są to Data Breakpoint, które nawiasem mówiąc pomogły mi ostatnio znaleźć pewien błąd który wydawałoby się brał się  znikąd.

Tak więc czym są owe Data Breakpoints? W przeciwieństwie do swoich „standardowych” kolegów nie zatrzymują wykonywania programu po osiągnięciu określonego miejsca, lecz w przypadku kiedy określony obszar pamięci ulega zmianie. Jak to działa? Posłużmy się prostym programem:

#include <iostream>

int g_Fred = 0;

void DoSomeBoringThings();

int main()
{
	std::cout << "Fred should be equal to Wilma\n";
	for(int Wilma = 0; Wilma < 10; Wilma++)
	{
		DoSomeBoringThings();
		std::cout << "Wilma equals " << Wilma << " and Fred equals " << g_Fred++ << "\n";
	}
	std::cin.get();
}

void DoSomeBoringThings()
{
	//some boring code
	++g_Fred;
	//some code written by a bored programmer
}

Program jest bardzo prosty i nie robi nic konkretnego więc skupmy się tylko na tym co istotne dla naszych rozważań. W funkcji DoSomeBoringThings modyfikujemy zmienną globalną. Jednak w pętli programista oczekuje że g_Fred nie zmienia się nigdzie ‚za kulisami’. Oczywiście w tak krótkim programie widać jak na dłoni co się dzieje, ale w większych projektach może być trudno znaleźć miejsce gdzie jakaś zmienna jest niespodziewanie modyfikowana. Tu z pomocą przychodzi nam tytułowe narzędzie, pozwalające zatrzymać program przy okazji modyfikowania fragmentu pamięci (więc także zmiennej). Jak tego użyć? Stawiamy zwyczajny breakpoint na początku funkcji main i uruchamiamy debugowanie. Kiedy program się zatrzyma zaznaczamy naszą zmienną (g_Fred) i wybieramy Debug->New Breakpoint->New Data Breakpoint… Pojawi nam się okno dialogowe New Breakpoint

W pole Addres wpisujemy adres pamięci, możemy użyć też wyrażenia z operatorem & (tak jak zrobiliśmy to w tym przykładzie). Byte Count oznacza ile bajtów od danego adresu mamy śledzić (przykładowo jest to 4 gdyż g_Fred jest typu int). Po wznowieniu, program zatrzyma się i zostaniemy przeniesieni w miejsce modyfikacji zmiennej.

Jak widać jest to nieskomplikowane narzędzie, jednakże może okazać się przydatne.

P.S. Post z jednodniowym poślizgiem, ale lepsze to niż wrzucenie wypełniacza tylko po to aby zmieścić się w zakładanym terminie.