10 Klasyfikacja do K grup, Sieć probabilistyczna GLM

1

plik wglm.tex, 26 listopada 2005, 28.02.2006, 5.12.09

10 Klasyfikacja do K grup, Sie´ c probabilistyczna GLM

Jest to jeden z najprostszych modeli sztucznych sieci neuronowych. Mo˙ze slu˙zy´c do ob- licze´n zagadnie´n formuÃlowanych jako modele liniowe (Linear Models) oraz tzw. uog´olnione modele liniowe (Generalized Linear Models ¹). Pomimo swej prostoty, sie´c GLM potrafi oblicza´c zagadnienia regresji wielokrotnej, liniowe zagadnienia najmniejszych kwadrat´ow, oraz klasyfikacje dla kilku (c) grup. Obliczenia sa bardzo szybkie, istnieje tylko jedno_, minimum optymalizujace zadane kryterium bÃl_, edu._,

10.1 Architektura sieci GLM

Architektura sieci GLM zaimplementowanej w pakiecie Netlab jest pokazana na Ry- sunku 10.1.

GLM neural network d=4 c=2

x1→ x2→

x3→

x4→

→ y1

→ y2

f(Σ) f(Σ)

d inputs # hidden = # output

Oznaczenia

d – l. wej´s´c, na rysunku d = 4 x = [x₁, . . . , x_d]^T – wektor wej´sciowy c – liczba neuron´ow warstwy ukrytej, i zarazem liczba wyj´s´c, na rysunku c = 2

w_j = [w_j1, . . . , w_jd]^T, j = 1, . . . , c wagi, z jakimi sumuje j-ty neuron b_j – bias j-tego neuronu

a_j – aktywacja j-tego neuronu, aj = x^T wj+ bj

f (.) – funkcja aktywacji, okre´slona na aktywacjach a_j, mo˙ze to by´c

’linear’, ’logistic’, ’softmax’

y_j = f (a_j) – wynik na wyj´sciu j-tego neuronu, j = 1, . . . , c

PoÃlaczenie mi_, edzy neuronami war-_, stwy ukrytej (wyj´sciowej) jest urucha- miane tylko w przypadku funkcji ak- tywacji softmax

Rysunek 10.1: Architektura sieci GLM o d = 4 wej´sciach i c = 2 wyj´sciach. Sie´c GLM ma tylko jedna warstw_, e ukryt_, a o c neuronach, kt´ora to warstwa jest jednocze´snie warstw_, a_, wyj´sciowa {myGLM.eps}_,

Oznaczenia:

x = (x₁, . . . , x_d)^T — wektor danych (input vector), zmienne obja´sniajace_, y = (y₁, . . . , y_c)^T — wektor wynik´ow (output vector), zmienne prognozowane

aj =^Pd_i=1wjixi+ bj, j = 1, . . . , c — aktywacje wej´sciowe neuron´ow warstwy ukrytej.

Wyj´scie yjjest obliczane za pomoca funkcji aktywacji f (.) wybieranej przez u˙zytkownika._, Argumentem funkcji f sa aktywacje a_{, j} (wynik sumowania, na rysznku oznaczony jako Σ), a rezultatem jest wynik y_j = f (a_j) stanowiacy odpowied´z sieci na sygnaÃl wej´sciowy x._,

1sa to takie modele, kt´ore przez prost_, a transformacj_, e daj_, a si_, e sprowadzi´c do modeli liniowych, por. [3]_,

Sie´c GLM w Netlabie ma wbudowane nastepuj_, ace trzy funkcje aktywacji:_,

• funkcja liniowa y_j = a_j

• funkcja logistyczna

y_j = 1

1 + exp( − a)

• funkcja softmax

exp(a_j)

P

j⁰exp(a_j⁰)

Ka˙zda z wymienionych funkcji aktywacji jest wÃla´sciwa dla innego zagadnienia:

funkcja liniowa – dla zagadnie´n regresyjnych typu y = b₀+ b₁x₁+ . . . + b_dx_d+ e

funkcja logistyczna – dla klasyfikacji do 2 grup; obliczona warto´s´c y ∈ [0, 1] wyra˙za prawdopodobie´nstwo, ˙ze wektor x nale˙zy do grupy okre´slanej jako ’klasa 1’

funkcja softmax – klasyfikacja danych do c grup (klas). Obliczone warto´sci wyj´sciowe (y₁, . . . , y_c) speÃlniaja warunki: 0 ≤ y_{, j} ≤ 1, ^Pc_j=1y_j = 1.

10.2 Implementacja sieci GLM w pakiecie Netlab

Istotna rol_, e odgrywaj_, a tu trzy funkcje: glm, glmfwd i glmtrain. S_, a one wywoÃlywane_, nastepuj_, aco:_,

net = glm(n in, n out, actfn);

Y = glmfwd(net, X); lub [Y, A] = glmfwd(net, X);

[net, options]= glmtrain(net, options, data, targets);

Funkcja GLM tworzaca struktur_, e net_, Pierwsza z funkcji tworzy obiekt ’net’.

Oczywi´scie zamiast nazwy ’net’ mo˙zna u˙zy´c dowolnej innej nazwy, np. ’my net’ lub ’jan’ : net = glm(n in, n out, actfn)

Znaczenie parametr´ow wej´sciowych:

n in : d, liczba neuron´ow warstwy wej´sciowej n out: c, liczba neuron´ow warstwy wej´sciowej

actfn : rodzaj funcji aktywacji: ’linear’, ’logistic’, ’softmax’.

Moga jeszcze wyst_, api´c dalsze parametry, np. prior, parametr opcyjny, wykorzystywany_, przy modelach Bayesowskich, okre´slajacy rozrzut inicjowanych wag._,

Utworzona struktura net ma nastepuj_, ace pola:_, type ’glm’

n in d, liczba wej´s´c (liczba cech danych) n out c, liczba wyj´s´c sieci

nwts liczba wszystkich wsp´oÃlczynnik´ow wagowych i bias´ow actfn funkcja aktywacji: string ’linear’, ’logistic’, ’softmax’

w1 tablica wag o wymiarze d × c b1 tablica bias´ow o wymiarze 1 × c

10.2 Implementacja sieci GLM w pakiecie Netlab 3

W przypadku u˙zywania zÃlo˙zonego modelu Bayesowskiego, struktura glmnet przewiduje dalsze pola zawierajace parametry Bayesowskie: prior i beta. Jednak my nie b_, edziemy_, zajmowa´c sie tym przypadkiem._,

Struktura ’net’ zainicjowana jak wy˙zej zawiera przyjete losowo warto´sci wag oraz bias´ow._, Np. Warto´sci wag zostaÃly zainicjowane jako randn(d, c)/sqrt(d − 1)

Funkcja GLMFWD sÃlu˙zaca pracy w trybie odtworzeniowym_,

Funkcja glmfwd pozwala utworzonej sieci pracowa´c w trybie odtworzeniowym, tzn. dla dostarczonych danych wej´sciowym dostarcza odpowiednie wyniki sieci.

WywoÃlywanie funkcji:

[Y, A] = glmfwd(net, X); lub Y = glmfwd(net, X);

Parametr X oznacza tu tablice danych X o wymiarze N ×d dla kt´orych chcemy otrzyma´c_, wyniki. Tablica ta zawiera N wektor´ow wej´sciowych:

X =







x1T

...

x_N^T





=







x11 x12 . . . x1d

... ... ... ...

x_{N 1} x_{N 2} . . . x_{N d}





 (10.1)

Jako rezultat dziaÃlania funkcji otrzymujemy tablice Y o wymiarze N × c oraz tablic_, e_, aktywacji A o wymiarze N × c. Sa one obliczane wedÃlug nast_, epuj_, acego wzoru:_,

|{z}A

N ×c

= X_|{z}

N ×d

|{z}W

d×c

+ 1_N ∗ b_|{z}

| {z1×c}

N ×c

, _|{z}Y

N ×c

= f (A)

| {z }

N ×c

,

gdzie W = [w₁, . . . , w_c] oraz b = [b₁, . . . , b_c] oznaczaja wagi i biasy zapami_, etane w polach_, net.w1 i net.b1 struktury ’net’; natomiast f (.) oznacza zadeklarowana funkcj_, e aktywacji._,

Oczywi´scie, je´sli wagi sa przypadkowe, to otrzymane wyniki b_, ed_, a te˙z przypadkowe._, Dlatego te˙z utworzona sie´c powinna by´c najpierw ’wyuczona’, inaczej m´owiac, wytreno-_, wana, za pomoca specjalnej procedury. Stanowi j_, a funkcja glmtrain._,

Funkcja glmtrain uczy sie´c rozpoznawania zadanych wzorc´ow.

Funkcja GLMTRAIN trenujaca sie´_, c wg. pr´obki uczacej_,

[net, options]= glmtrain(net, options, data, targets);

Parametry wej´sciowe:

options. Praca modelu, a przede wszystkim otrzymywane jako rezultat tej pracy wyniki, sa uwarunkowane tzw. opcjami. Opcje te s_, a zapami_, etywane w tablicy options(1:18)_, ². Po uruchomieniu systemu Netlab jest dostepna tablica foptions(1:18), kt´ora ma zainicjowane_, warto´sci domy´slne. U˙zytkownik mo˙ze skopiowa´c te tablic_, e do swojej tablicy options, a_, nastepnie – ju˙z na gruncie swoich ’options’ – nada´c poszczeg´olnym elementom odpowiednie_, warto´sci.

2Najwa˙zniejsze znaczenia warunk´ow zakodowanych w tablicy options sa opisane w komentarzu do funkcji_, glmtrain

Najwa˙zniejsze znaczenia warunk´ow zakodowanych w tablicy options sa opisane w ko-_, mentarzu do funkcji glmtrain. Mo˙zna r´ownie˙z otrzyma´c opis tych warunk´ow przez rozkaz help foptions. Najcze´sciej wykorzystywane s_, a nast_, epuj_, ace opcje dotycz_, ace procesu uczenia_, (trenowania):

options(1) – warto´s´c 1 oznacza drukowanie bÃledu E po wykonaniu ka˙zdej iteracji; 0 ozna-_, cza pomijanie tych wydruk´ow,

options(2) – wymagana dokÃladno´s´c wag;

options(3) – wymagana dokÃladno´s´c bÃledu E;_,

options(8) – zwracana warto´s´c bÃledu (po wykonaniu oblicze´n),_, options(14) – maksymalna liczba iteracji.

data . Tablica ’data’ powinna zawiera´c tzw. pr´obke ucz_, ac_, a o postaci tablicy X_, N ×d (zdefi- niowanej wzorem [10.1]). Liczba N odpowiada w tym przypadku liczbie wektor´ow uczacych_, x_i wchodzacych w skÃlad pr´obki ucz_, acej._,

targets. Jest to tablica pokazujaca wzorce, kt´ore sie´c ma si_, e nauczy´c rozpoznawa´c. Tablica_, ta jest postaci podobnej jak tablica ’data’, z tym, ˙ze zawiera tylko c kolumn. Oznaczmy te tablic_, e symbolem T = (t_{, ij}).

T =







t1T

...

tNT





=







t11 t12 . . . t1c

... ... ... ...

tN 1 tN 2 . . . tN c





 (10.2)

W zagadnieniach regresyjnych tablica T zawiera warto´sci numeryczne, natomiast w za- gadnieniach klasyfikacyjnych – zmienne zero-jedynkowe, okre´slajace przynale˙zno´s´c do posz-_, czeg´olnych klas.

Pr´obka uczaca skÃlada si_, e z dwu tablic: tablicy ’data’, zawieraj_, acej zmienne obja´sniaj_, ace,_, oraz tablicy ’targets’ zawierajacej warto´sci docelowe. Jest odpowiednio´s´c mi_, edzy wierszami_, obu tablic. Ka˙zdy wiersz t_j^T tablicy ’targets’ stanowi zbi´or wynik´ow, jakich sie´c ma do- starczy´c jako odpowied´z na wektor zmiennych obja´sniajacych x_{, j}^T.

data =







x₁^T ...

xNT





 targets =







t₁^T ...

tNT







rozmiar N × d rozmiar N × c

Proces trenowania. Podczas uczenia sie´c poprawia iteracyjnie swoje wagi w ten spos´ob, aby wyniki dostarczone przez sie´c (tj. tablica Y) mo˙zliwie maÃlo r´o˙zniÃly sie od_, wektora warto´sci po˙zadanych (docelowych) danych w tablicy T (’targets’). W tym celu_, nale˙zy okre´sli´c bÃlad E oraz algorytm aktualizuj_, acy wagi. Zagadnienia okre´slania bÃl_, edu_, sieci i jego minimizacji sa omawiane w rozdziale 3._,

Skuteczno´s´c wytrenowania sieci sprawdzamy na oddzielnej pr´obce nazywanej pr´obka_, testowa, dla kt´orej oblicza si_, e wielko´s´c bÃl_, edu. Dla zagadnie´n klasyfikacyjnych wielko´s´c_, bÃledu sprawdza si_, e na tzw. macierzy pomieszania (ang. confusion matrix, na podstawie_, kt´orej liczy sie liczb_, e element´ow sklasyfikowanych poprawnie._,

10.3 Obliczanie bÃledu sieci_, 5

10.3 Obliczanie bÃl edu sieci

[e, y, a]= glmerr(net, x, t);

Wyniki ’y’ i ’a’ sa opcyjne. Mo˙zna wywoÃlywa´c tylko: e= glmerr(net, x, t);_, Parametry wej´sciowe:

net – struktura, utworzona funkcja ’glm’ i ewtl. zmodyfikowana przez ’glmtrain’,_, x – tablica danych, wymiaru N × d,

t – tablica warto´sci docelowych, wymiaru N × c.

Wyniki:

e – tablica wymiaru N × 1, bÃlad sieci, zale˙zy od przyj_, etej funkcji aktywacji_, opcyjnie:

y – wyniki sieci podawane na wyj´sciu, czyli tablica Y_{N ×c} a – aktywacje neuron´ow warstwy ukrytej, tablica

AN ×c= [ ˜X ˜w1, . . . , ˜X ˜wc] = X ˜W,

znak tildy oznacza tu odpowiednio poszerzone tablice X – o kolumne jedynek, oraz W –_, o biasy neuron´ow:

X = [X, 1˜ N], W = [ ˜˜ w1, . . . , ˜wc] =

"

w1 . . . wc

b₁ . . . b_c

#

.

BÃlad e mo˙ze by´c zapami_, etany w tablicy options jako element options(8)._, A oto tre´s´c funkcji obliczajacej bÃl_, ad sieci:_,

[y, a] = glmfwd(net, x);

switch net.outfn

case ’linear’ % Linear outputs edata = 0.5*sum(sum((y - t).^2));

case ’logistic’ % Logistic outputs

edata = - sum(sum(t.*log(y) + (1 - t).*log(1 - y)));

case ’softmax’ % Softmax outputs edata = - sum(sum(t.*log(y)));

otherwise

error([’Unknown activation function ’, net.outfn]);

end

10.4 Obliczanie wag neuron´ ow

Algorytm obliczania wag zale˙zy od przyjetej funkcji aktywacji. W Netlabie dopuszcza_, sie nast_, epuj_, ace funkcje aktywacji: liniowa, logistyczna, softmax._,

Wprowadzamy tu specyficzne oznaczenie wag neuron´ow: bedziemy je oznacza´c sym-_, bolem β (a nie symbolem w, jak to robimy w innych rozdziaÃlach).

Dla c neuron´ow warstwy ukrytej szukamy c wektor´ow wagowych β_j, gdzie β_j = [β₁, . . . , β_d, β₀]^T, j = 1, . . . , c.

Ostatni, (d + 1)-szy element oznacza bias j-tego neuronu.

Liniowa funkcja aktywacji

W przypadku liniowej funkcji aktywacji funkcje bÃl_, edu E wyznacza si_, e na zasadzie naj-_, mniejszych kwadrat´ow r´o˙znic (metoda Least Squares). Aby wyznaczy´c najlepsze wagi i biasy, jest jednorazowo rozwiazywany ukÃlad r´owna´n liniowych_,

XB = T (10.3)

gdzie: X_{N ×(d+1)} – tablica danych z dodatkowa kolumn_, a ’jedynek’ na ostatnim miejscu,_, B_(d+1)×c – wagi neuron´ow, kolumnami, ostatnia skÃladowa kolumny zawiera bias danego neuronu:

B = [β₁, . . . , β_c] =

"

w₁ . . . w_c b1 . . . bc

#

T_{N ×c} – tablica warto´sci docelowych.

W przypadku sieci GLM mamy tyle neuron´ow co wyj´s´c, stad H ≡ c._,

Wypisany wy˙zej ukÃlad r´owna´n (10.3) jest nadokre´slonym ukÃladem r´owna´n liniowych na wagi i biasy sieci. UkÃlad ten rozwiazuje si_, e w Netlabie przez operacj_, e matrix left divide._,

Logistyczna funkcja aktywacji

Tutaj jako funkcje bÃl_, edu przyjmujemy E = − logL; patrz rozdziaÃl 3 notatek i skrypt_, w paragrafie 4.3 powy˙zej. Zale˙zno´sci miedzy danymi wej´sciowymi (X), wagami sieci (B) i_, warto´sciami docelowymi (T) sa nieliniowe, ale – przy logistycznej funkcji aktywacji – pod-_, padaja pod tzw. uog´olniony model liniowy (GLM) opisany dalej w sekcjach 4.7 i 4.8. Al-_, gorytm znajdywania wag w takim modelu sprowadza sie do sekwencyjnego rozwi_, azywania_, ukÃladu r´owna´n liniowych ”wa˙zonych” (wagi indywidualne dla ka˙zdego osobnika) (10.4), w k´orych jako niewiadoma wystepuje wektor wag β_{, j} charakteryzujacy j-ty neuron_, ³, j = 1, . . . , c.

Rozwa˙zymy najpierw przypadek, gdy liczba klas = 2, a c = 1. Jest to sto- sunkowo czesto stosowany przypadek, gdy˙z przy liczbie klas > 2 stosujemy raczej funkcj_, e_, softmax.

Przy dw´och klasach mamy: B_(d+1)×1 = β_(d+1)×1, czyli sie´c skÃlada sie z jednego neuronu_, i mamy tylko jeden wektor wag β = (β₁, . . . , β₀)^T; r´ownie˙z T = t_{N ×1}.

3mamy tu dwa rodzaje wag: osobnicze, okre´slajace ’wag_, e’ ka˙zdego wiersza tablicy X przy wyznaczaniu_, warto´sci docelowych , T, oraz wagi neuron´ow u˙zywanych do wyznaczania aktywacji neuronu w odpowiedzi na sygnaÃl wej´sciowy xⁿ.

10.5 ModuÃly demonstracyjne do modelu GLM 7

Niech Π_{N ×1} = (π₁, . . . , π_N)^T oznacza wektor prawdopodobie´nstw, ˙ze i-ty osobnik nale˙zy do klasy ’1’ (’sukces’ w rozkÃladzie binomialnym, przy osobnikach zr´o˙znicowanych ze wzgledu na zmienne obja´sniaj_, ace). Prawdopodobie´nstwa te b_, edziemy szacowa´c ite-_, racyjnie. Najpierw przyjmiemy jaka´s warto´s´c pocz_, atkow_, a p, a nast_, epnie, w kolejnych_, iteracjach, warto´s´c te b_, edziemy ’poprawia´c’ (m´owi_, ac dokÃladniej, warto´s´c ’poprawion_, a p_, bedziemy otrzymywa´c – z procedury glmerr)._,

Warto´s´c Π posÃlu˙zy do skonstruowania ukÃladu r´owna´n (rLS)

W^1/2Xβ^(k+1)= W^1/2z, (10.4)

gdzie: W= diag(w11, . . . , wN N) oznacza wagi i-tego osobnika (r´ownania), przy czym w_ii= [p_i(1 − p_i)]^1/2, i = 1, . . . , N ,

z_{N ×1}= log(p./(1 − p)) + (t − p)./[p(1 − p)],

wektor pomocniczy obliczony z szacunku p w (k)-tej iteracji.

Trenowanie sieci odbywa sie w postaci wsadowej (batch) na podstawie caÃlej pr´obki_, uczacej, tj. wszystkich element´ow {x_{, i}, t_i}, (i = 1, . . . , N) – przy u˙zyciu wzoru (10.4).

Inicjacja (krok k = 0) polega na wstepnym przyj_, eciu p = (t + 0.5)/2. Z warto´sci tych_, wyznaczamy wstepne warto´sci W i z; nast_, epnie rozwi_, azujemy ukÃlad r´owna´n ze wzgl_, edu na_, wektor β. Znajomo´s´c β (wag sieci) pozwala nam zaktualizowa´c warto´s´c p.

Czynno´sci te powtarzamy dla k = 1, . . . , a˙z osiagniemy ˙z_, adan_, a dokÃladno´s´c, lub wy-_, czerpiemy zadeklarowana liczb_, e iteracji._,

W ka˙zdej iteracji: (a) Przy znanym p wyznaczamy W oraz z.

(b) Majac te warto´sci, rozwi_, azujemy ukÃlad r´owna´n 10.4 – robimy to przez operacj_, e_, matrix left divide; nastepnie wyznaczamy bÃl_, ad sieci i nowe warto´sci p._,

Przypadek, gdy liczba klas c jest wieksza od 2. Rozbijamy zagadnienie na c_, niezale˙znych problem´ow z dwoma klasami i rozwiazujemy ka˙zdy problem oddzielnie, jak w_, przypadku 2 klas.

Warto´sci stopu

Funkcja aktywacji softmax

W przypadku aktywacji softmax i options(5)=1, to powtarzamy c-krotnie (tyle razy, ile jest wyj´s´c) postepowanie stosowane przy aktywacji logistycznej. Ostateczne wyniki (p-stwa_,

´a posteriori) sa skalowane tak, aby sumowaÃly si_, e do jedno´sci. W przypadku options(5)=0_, stosujemy dokÃladne obliczenia Hesjanu.

10.5 ModuÃly demonstracyjne do modelu GLM

Pakiet Netlab oferuje moduÃly demglm1 i demglm2 ilustrujace klasyfikacj_, e do dw´och i trzech_, grup przy u˙zyciu sieci GLM korzystajacej z funkcji aktywacji ’logistic’ (2 grupy danych) i_,

’softmax’ (3 grupy danych).

10.5.1 ModuÃl demglm1

W module demglm1 generuje sie najpierw dwie grupy danych dwuwymiarowych. Wy-_, generowne dany zostaja wykre´slone na pÃlaszczy´znie, ka˙zda grupa innym kolorem. Dane_,

zostaja podzielone na pr´obk_, e ucz_, ac_, a i pr´obk_, e testow_, a czyli sprawdzaj_, ac_, a. Sie´c glm ma_, tutaj architekture_,

input: d=2 → hidden: H=1 → output: c=1

Macierz warto´sci docelowych skÃlada sie tylko z jednej kolumny t_{, N ×1}, zawierajacej ’1’,_, gdy wygenerowany wektor x^T nale˙zy do pierwszej grupy, oraz warto´s´c ’0’, gdy nale˙zy do grupy drugiej.

Nastepnie tworzymy struktur_, e ’net’ za pomoc_, a konstruktora glm (deklaruj_, ac funkcj_, e_,

’logistic’ jako funkcje aktywacji) i trenujemy utworzon_, a sie´c za pomoc_, a funkcji glmtrain._, Jako rezultat trenowania otrzymujemy wagi sieci zaadaptowane do rozpoznawania klas danych. Samo rozpoznawanie jest wykonywane za pomoca funkcji glmfwd. Dla danego_, wektora x, kt´ory ma by´c rozpoznany, sie´c oblicza najpierw aktywacje (znak tildy oznacza_, odpowiednio poszerzone wektory x- poszerzony o warto´s´c 1, i w - poszerzony o bias)

a = a(x) = ˜x^Tw,˜ a nastepnie sw´oj wynik_,

y = y(x) = 1

1 + exp {−a(x)}.

Na tym samym wykresie, na kt´orym wyrysowane punkty indywidualne wygenerowanej pr´obki, mo˙zna r´ownie˙z wykre´sli´c powierzchnie funkcji y(t|x, w) i jej warstwice. W szcze-_, g´olno´sci interesujace s_, a warstwice dla y = 0.1, 0.5, 0.9. Odpowiedni fragment skryptu_, wykre´slajacy te warstwice przedstawiono pony˙zej._,

% fragment skryptu glmdem1.m,

% wykresla kontury i powierzchnie ’mesh’ i ’surf’

x = -4.0:0.2:5.0; y = -4.0:0.2:5.0; % jak g"esta ma byc siatka

[X, Y] = meshgrid(x,y); % utworzone X, Y sa dwuwymiarowe X = X(:); Y = Y(:); % teraz X, Y s"a kolumnami

Z = glmfwd(net,[X Y]); % [X Y] to tablica danych o 2 kolumnach Z = reshape(Z, length(x), length(y));

v = [0.1 0.5 0.9]; % wysoko"sci na ktorych przekroje [c, h] = contour(x, y, Z, v); % wykreslanie konturow na istniejacym rysunku title(’Model GLM’, ’Fontsize’,12)

set(h, ’linewidth’, 3)

clabel(c,h,’fontsize’,15,’fontweight’,’bold’,’color’,’r’,...

’labelspacing’,72,’rotation’,0) % opisywanie konturow fh4=figure, surf(x,y,Z) % Wykreslanie powierzchni pokrytej plytkami fh5=figure, mesh(x,y,Z) % Wykreslanie powierzchni zaznaczonej liniami

Jako´s´c klasyfikacji mo˙zna sprawdzi´c za pomoca tzw. macierzy pomieszania (confusion_, matrix) – otrzymamy wtedy tablice 2×2 pokazuj_, ac_, a dla i=0,1 ile element´ow klasy i zostaÃlo_, rozpoznanych jako nale˙zacych do klasy 2 (oznaczanej jako ’0’), a ile – do klasy 1._,

Taka tablic_, e mo˙zemy otrzyma´c za pomoc_, a funkcji conffig (wywoÃlanie: conffig(y test,t test);)._,

10.5.2 ModuÃl demglm2

W modelu tym por´ownuje sie wyniki otrzymane z sieci neuronowej z wynikami otrzy-_, manymi Bayesowska metod_, a prawdopodobie´nstw a posteriori. Jest to przedstawione za_, pomoca interesuj_, acej grafiki._,

10.6 Uczenie sieci: liniowa funkcji aktywacji 9

10.6 Uczenie sieci: liniowa funkcji aktywacji

Liniowa funkcja aktywacji prowadzi do funkcji bÃledu, kt´ora jest kwadratow_, a funkcj_, a wag. Szu-_, kamy wtedy wag , kt´ore daja minimum nast_, epuj_, acej formy kwadratowej {normalEq}:_,

W₀: min

W (Y − XW)^T(Y − XW), (10.5)

gdzie W_(d+1)×c – wagi i biasy podstawione do jednej tablicy, Y_{N ×c} – tablica zawierajaca wyniki sieci_,

X_{N ×(d+1)} – tablica danych ’data’ poszerzona o kolumne jedynek I_{, N}.

Warunek (10.5) prowadzi do tzw. ukÃladu r´owna´n normalnych, bed_, acych ukÃladem liniowych_, r´owna´n na elementy macierzy W, dajacych rozwi_, azanie w jednym kroku._,

Alternatywnie mo˙zna otrzyma´c rozwiazanie przez pseudoodwrotno´s´c. Rozpatrujemy wtedy_, ukÃlad r´owna´n

XW = Y,

kt´ore rozwiazujemy przez pseudoodwrotno´s´c (funkcja pinverse w Matlabie)._,

Niech X^† oznacza pseudoodwrotno´s´c macierzy X. Wtedy W = X^†Y daje nam szukane rozwiazanie._,

10.7 Uczenie sieci: logistyczna funkcja aktywacji

Funkcje bÃledu dla zagadnie´n klasyfikacji do 2 lub c > 2 klas zostaÃly przedstawione w_, rozdziale 3. Sa to wzory wynikaj_, aj_, ace z odpowiednich modeli probabilistycznych (bino-_, mialnych) i zwiazan_, a z tymi modelami funkcj_, a wiarogodno´sci obserwacji. Funkcja bÃl_, edu_, E oblicza sie wtedy jako logarytm z wiarogodno´sci L wzi_, ety ze znakiem przeciwnym. Wy-_, nikajace st_, ad wzory na bÃl_, edy (entropia krzy˙zowa) zale˙z_, a od wag w, kt´ore nale˙zy wyznaczy´c_, iteracyjnie w procesie uczenia (treningu) sztucznej sieci neuronowej.

Je˙zeli funkcja bÃledu E jest r´o˙zniczkowalna, to w celu minimizacji bÃl_, edu mo˙zemy sto-_, sowa´c metody gradientowe. Zasada jest wtedy taka, ˙ze dysponujac wektorem w_, ^(k)otrzyma- nym w k-tym kroku, “poprawiamy” ten wektor poruszajac si_, e o maÃl_, a odlegÃlo´s´c w kierunku,_, w kt´orym funkcja E maleje najbardziej, czyli w kierunku ujemnego gradientu.

W ten spos´ob w kolejnych krokach {k}, k = 0, 1, . . . mamy nadzieje zbli˙za´c si_, e do_, minimum funkcji E.

Na poczatku algorytmu musimy poda´c jakie´s wst_, epne przybli˙zenie w_, ⁽⁰⁾. Czesto jest to_, przybli˙zenie wygenerowane losowo.

W dalszych krokach⁴ ciag kolejnych przybli˙ze´n jest konstruowany wg zasady: maj_, ac_, wektor w^(k), konstruujemy w^(k+1) wedlug wzoru

w^(k+1) = w^(k)− η∂E

∂w

¯¯

¯¯

w=w^(k).

Wsp´oÃlczynnik η nosi nazwe wsp´oÃlczynnika uczenia (learning rate)._,

Jednak nie za bardzo wiadomo, jaka ma by´c warto´s´c η. Czesto warto´s´c t_, e wyznacza si_, e_, eksperymentalnie.

Jedna z uznawanych metod numerycznych, gwarantuj_, acych zbie˙zno´s´c (chocia˙z mo˙ze to_, by´c tylko zbie˙zno´s´c lokalna) jest metoda Newtona - Raphsona. Metoda ta poleca konstru- owa´c (k + 1)sze przybli˙zenie wag wg wzoru:

w^(k+1)= w^(k)− H⁻¹∂E

∂w.

4przedstawimy dalej algorytm opisany w ksia˙zce Nabney’a [4]_,

Metoda ta wymaga odwracania hesjanu, co przy wiekszych rozmiarach wag w mo˙ze by´c_, kÃlopotliwe, r´ownie˙z dlatego, ˙ze wymagana jest dodatnia okre´slono´s´c macierzy H.

Fisher zaproponowaÃl metode, kt´ora zamiast hesjanu H podstawia warto´s´c oczekiwan_, a_, tej macierzy (metoda ta jest nazywana Fisher’s scoring method⁵, u nas E = −logL, gdzie L jest wiarygodno´scia pr´oby, a wektor w peÃlni rol_, e wektora parametr´ow ) {wz´or FisherScor}_,

w^(k+1) = w^(k)− {E[H]}⁻¹∂E

∂w. (10.6)

Je´sli Funkcja aktywacji f jest funkcja liniow_, a, to Hesjan dla tak okre´slonego modelu_, liniowego wynosi X^TX [4], str 129–132. Hesjan dla przypadku f=’logistic’ jest pokazany ni˙zej.

Rozpatrzmy zagadnienie klasyfikacyjne dla dw´och klas okre´slonych binarnie, gdy dla danego wektora danych x sie´c ma poda´c wynik yⁿ(x) okre´slajacy przynale˙zno´s´c do_, grupy (klasy) 1 (0 ≤ yⁿ ≤ 1. Prawdziwa przynale˙zno´s´c do klasy jest okre´slona podana_, warto´scia t_, ⁿ∈ {0, 1}:

xⁿ → yⁿ ma przybli˙za´c tⁿ ZakÃladamy, ˙ze funkcja aktywacji f jest funkcja logistic_,

Niech Π oznacza macierza diagonaln_, a: Π = diag{y_, ⁿ(1 − y)ⁿ}.

Mamy (wzory cytowane za [4], str 129–132)

• Gradient ∂E/∂w = X^TΠe, gdzie eⁿ = (yⁿ− tⁿ)/f⁰(aⁿ), f - funkcja aktywacji, aⁿ – aktywacja neuronu w odpowiedzi na sygnaÃl wej´sciowy xⁿ,

• Hesjan dla tego samego modelu r´owna sie X_, ^TΠX,

Wypisane wy˙zej wzory pozwalaja na prosty algorytm oblicze´n, zastosowany w pakiecie_, Netlab. Dla posiadanej warto´sci w^(k) tworzymy zmienna_,

z_k= Xw^(k)+ e.

Nastepnie rozwi_, azujemy iteracyjnie wa˙zone liniowe zagadnienia najmniejszych kwadrat´ow_, (iterrative reweighted least squares equations) czyli nastepuj_, acy ukÃladu r´owna´n normal-_, nych (wektor w jest tu wektorem kolumnowym rozmiaru (d + 1) × 1, a X macierza danych_, poszerzona o kolumn_, e jedynek): {wz´or iksPi}_,

(X^TΠkX)w^(k+1)= X^TΠkzk. (10.7)

UkÃlad wypisany powy˙zej jest r´ownowa˙zny zwykÃlemu ukÃladowi r´owna´n normalnych, je´sli podstawimy za macierz danych tablice X_, ^TΠ^1/2_k , a za zmienna zale˙zn_, a tablic_, e Π_, ^1/2_k z_k.

Warto´sci Πk i zk zmieniaja si_, e w ka˙zdej iteracji, co zaunaczono wska´znikiem ’k’._, W przypadku klasyfikacji do wiecej grup i u˙zywania funkcji aktywacji ’softmax’_, metoda minimalizacji bÃledu zale˙zy od options(5); je˙zeli options(5)=1, to stosuje si_, e metod_, e_, IRLS jak w przypadku f=’logistic’ (Nabney [4] podaje r´ownie˙z wzory na single output so- ftmax model. Nie sa one du˙zo bardziej skomplikowane ani˙zeli te, podane wy˙zej dla modelu_, logistycznego), w przeciwnym przypadku jest stosowana dokÃladna metoda numeryczna Newtona-Raphsona z obliczaniem odwrotno´sci hesjanu.

W wyniku trenowania otrzymuje sie aktualizacj_, e wag i bias´ow (og´olnie: wektora wag w)_, oraz nowa warto´s´c bÃl_, edu E. Ten ostatni jest zapami_, etany jako options(8)._,

5Fisher w swoim podej´sciu rozwa˙zal logarytm wiarygodno´sci l = logL przedstawiany jako funcja para- metru θ i szukaÃl maksimum funkcji l ze wzgledu na parametr θ. Funcja l(θ) wyznaczona dla N-elementowej_, pr´oby o niezale ˙nych elementach dawaÃla sie przedstawi´c jako suma_, P_N

i=1li(θ))

10.8 Metoda GLM – uog´olniony model liniowy 11

10.8 Metoda GLM – uog´ olniony model liniowy

W tej sekcji zostanie przedstawiona teoria uog´olnionych modeli liniowych GLM rozwinieta_, na gruncie statystyki matematycznej. Model ten (tj. GLM) obejmuje jako szczeg´olny przypadek zagadnienie klasyfikacji do 2 grup danych, gdy p-stwo klasyfikacji π zale˙zy od obserwowanych dla ka˙zdego osobnika zmiennych obja´sniajacych X_, 1, . . . , Xd.

10.8.1 Rodzina wykÃladnicza rozkÃlad´ow probabilistycznych Definicja rodziny wykÃladniczej

RozkÃlad f (y; θ) nale˙zy do rodziny wykÃladniczej rozkÃlad´ow, je´sli daje sie zapisa´c w postaci:_, {wykladniczy}

f (y; θ) = exp [ a(y)b(θ) + c(θ) + d(y)], (10.8) gdzie funkcje b,c,d sa r´o˙zniczkowalne conajmniej dwa ray._,

Je˙zeli a(y) ≡ y, to m´owimy, ˙ze b(θ) jest parametrem naturalnym rozkÃladu, a posta´c (10.8) jest postacia kanoniczn_, a rozkÃladu._,

Zaznaczmy tutaj, ˙ze wiele znanych i u˙zywanych w praktyce rozkÃlad´ow probabilistycz- nych nale˙zy do tej rodziny; miedzy nimi s_, a r´ownie˙z rozkÃlad binarny (zero-jedynkowy, Ber-_, noulliego), rozkÃlad dwumianowy, Poissona, multinomialny, gamma, i rozkÃlad normalny.

Przy rozpatrywaniu zagadnie´n klasyfikacyjnych interesuje nas najbardziej rozkÃlad zero- jedynkowy (Bernoulliego), bed_, acy szczeg´olnym przypadkiem rozkÃladu dwumianowego (dla_, n = 1).

PrzykÃlad: RozkÃlad dwumianowy Y ∼ binomial(n, π). Tutaj π, prawdopodobie´nstwo

’sukcesu’, jest interesujacym nas paramerem, a n, dÃlugo´s´c serii, jest dane. Dla n = 1_, otrzymujemy rozkÃlad binarny.

Funkcja rozkÃladu p-stwa daje sie zapisa´c w postaci kanonicznej jako (y = 0, 1, . . . , n):_, f (y; π) = exp^hylogπ − ylog(1 − π) + n(1 − π) + log

Ãn y

!i

Tak wiec mamy tu: b(π) = logπ − log(1 − π) = log[π/(1 − π)]._,

Jest to zarazem parametr naturalny dla rozkÃlad´ow dwumianowego i binarnego (zero–

jedynkowego).

WÃlasno´sci rozkÃlad´ow rodziny wykÃladniczej

E[a(Y )] = −c⁰(θ)/b⁰(θ). (10.9)

V ar[a(Y )] = b⁰⁰(θ)c⁰(θ) − c⁰⁰(θ)b⁰(θ)

[b⁰(θ)]³ . (10.10)

Warto´sc oczekiwana i wariancja sa tu liczone ze wzgl_, edu na rozkÃlad zmiennej losowej Y ._, PrzykÃlad: RozkÃlad dwumianowy c.d. Dla rozkÃladu dwumianowego mamy:

θ = π, a(Y ) ≡ Y, b(π) = log[π/(1 − π)], c(π) = n log(1 − π) b⁰(θ) = 1

π(1 − π), b⁰⁰(θ) = π − (1 − π)

π²(1 − π)² , c⁰(θ) = −n

1 − π, c⁰⁰(θ) = −n (1 − π)². Skad – poniewa˙z a(Y ) ≡ Y – znane wzory: E(Y ) = nπ, V ar(Y ) = nπ(1 − π)._,

Wiarygodno´s´c l(θ; y) = log L i jej pochodne dla rodziny wykÃladniczej

Niech l(θ; y) = log L oznacza logarytm funkcji wiarogodno´sci. Dla rodziny wykÃladniczej mamy

l(θ; y) = a(y)b(θ) + c(θ) + d(y)

Aby wyznaczy´c parametr dajacy maximum wiarygodno´sci, nale˙zy obliczy´c pochodn_, a wzgl_, edem_, θ i przyr´owna´c ja do zera (akÃladamy, ˙ze funkcja ta jest dostatecznie regularna, a maximum_, nie le˙zy na brzegu dziedziny funkcji pochodnej). Dla rodziny wykÃladniczej definiujemy:

U(θ; y) = ∂l(θ; y)

∂θ = a(y)b⁰(θ) + c⁰(θ).

Definicja Funkcja U nosi nazwe score statistics._, Warto´s´c oczekiwana funkcji U : E(U)= 0. Dow´od:

E(U) = {E[a(Y )]}b⁰(θ) + c⁰(θ) = [− c⁰(θ)/b⁰(θ)] b⁰(θ) + c⁰(θ) = 0.

Wariancja funkcji U : V ar(U) = E(U²) = − E(U⁰), gdzie U⁰ = ∂U/∂θ.

Wariancje U mo˙zemy r´ownie˙z obliczy´c z og´olnego wzoru jako wariancj_, e transformowanej_, liniowo zmiennej a(y):

I = V ar(U) = [b⁰(θ)²] V ar[a(Y )] = b⁰⁰(θ)c⁰(θ)

b⁰(θ) − c⁰⁰(θ)

Zauwa˙zmy, ˙ze I zale˙zy jedynie od parametr´ow przyjetego modelu, nie zale˙zy natomiast od_, warto´sci oczekiwanych Y .

Definicja. Wariancja U nosi nazwe informacji._,

Jako pochodna U⁰ = ∂U/∂θ (dokÃladnie: V ar(U) = − E(U⁰)) wyznaczona w punkcie ˆθ dostarcza ona informacji, jak szybko zmienia sie funkcja wiarogodno´sci (a wÃla´sciwie jej_, logarytm) w okolicy maksimum. Je´sli pochodna ta zmienia sie bardzo powoli (funkcja_, jest pÃlaska), to precyzja wyznaczonego estymatora jest maÃla, a wariancja wyznaczonego estymatora du˙za.

PrzykÃlad: RozkÃlad dwumianowy c.d. Dla rozkÃladu dwumianowego mamy:

U = a(Y )b⁰(θ) + c⁰(θ) = y

π(1 − π)− −n 1 − π.

Biorac pod uwag_, e, ˙ze dla r. dwumianowego V ar(Y ) = nπ(1 − π), otrzymujemy_, V ar(U) = nπ(1 − π)

π²(1 − π)² = n

π(1 − π) = I . 10.8.2 Og´olna koncepcja GLM

Metoda GLM jest oparta na statystycznej teorii uog´olnionych modeli liniowych wprowa- dzonych przez Neldera i Wedderburna w 1972 roku (por [7]). Teoria ta oraz referencje do innych monografii i prac na temat GLM sa podane m.in. w ksi_, a˙zce Annette Dobson [3]._,

Metoda GLM uog´olnia statystyczna teori_, e modelu regresyjnego (czyli wyznaczania_, zale˙zno´sci (regresji) zmiennej Y od zmiennych obja´sniajacych X_{, 1}, . . . , X_d).

U podstaw uog´olnionego modelu liniowego le˙zy zaÃlo˙zenie, ˙ze obserwowane zmienne losowe Y1, . . . , YN zale˙za od indywidualnych parametr´ow θ_, 1, . . . , θN, kt´ore to parametry sa_, w praktyce nieznane. Co gorsza – parametry te sa niemo˙zliwe do wyestymowania, je˙zeli_,

10.8 Metoda GLM – uog´olniony model liniowy 13

dysponujemy tylko jedna obserwacj_, a dla ka˙zdego Y_{, i}. Ponadto, przy wiekszych pr´obkach,_, parametr´ow do wyestymowania byÃloby zbyt wiele.

W tej sytuacji powstaÃla teoria uog´olnionych modeli liniowych GLM, nazywana r´ownie˙z GLIM, od Generalzed Linear Models. Teoria ta stara sie_,

a) zmniejszy´c liczbe parametr´ow opisuj_, acych dan_, a zbiorowo´s´c_,

b) przy mo˙zliwie maÃlej liczbie parametr´ow wymodelowa´c zale˙zno´sci nieliniowe obserwo- wanej zmiennej Y od zmiennych obja´sniajacych X_{, 1}, . . . , X_d.

Metoda GLM stara sie scharakteryzowa´c rozkÃlady obserwowanych zmiennych losowych_, Y₁, . . . , Y_N mniejsza liczb_, a parametr´ow β = (β_{, 1}, . . . , β_d)^T, (d < N) takich, ˙ze kombinacja liniowa przyjmowanych parametr´ow β₁, . . . , β_d jest pewna funkcj_, a warto´sci oczekiwanej_, i-tej obserwowanej warto´sci zmiennej Y_i.

Innymi sÃlowy: Idea uog´olnionego modelu liniowego jest ”that the image of the mean_, response by a given link function can be modelled via a linear relationship.” [6].

Niech Y₁, . . . , Y_N bed_, a obserwowanymi warto´sciami pr´obkowymi oznaczaj_, acymi warto´sci_, zmiennej Y zaobserwowane dla N r´o˙znych osobnik´ow. ZakÃladamy, ˙ze:

• Realizacje Y₁, . . . , Y_N sa niezale˙zne._,

• RozkÃlad Y_i, (i = 1, . . . , N ) opisuje sie rozkÃladem f (y; θ_{, i}), tj. tym samym rozkÃladem dla ka˙zdego i, ale z parametrem θ_i kt´ory mo˙ze by´c indywidualny dla osobnika nr. i,

• RozkÃlad f (yi; θi) nale˙zy do rodziny wykÃladniczej rozkÃlad´ow i jest postaci kanonicznej, co oznacza, ˙ze daje sie zapisa´c w postaci: {wykladniczyCan}_,

f (y_i; θ_i) = exp [ y_ib(θ_i) + c(θ_i) + d(y)]. (10.11)

• RozkÃlad Yi zale˙zy od obserwowanych zmiennych obja´sniajacych X_, 1, . . . , Xd kt´orych warto´sci sa znane: tzn. dla ka˙zdego Y_{, i}znamy odpowiadajace warto´sci x_, ^T_i = (x_i1, . . . , x_id).

Zale˙zno´s´c ta jest postaci

g(µi) = x^T_i β. (10.12)

gdzie:

µ_i oznacza warto´s´c oczekiwana zmiennej losowej Y_{, i}: tj. mamy E(Y_i) = µ_i. g(.) jest tzw. funkcja–Ãl_, acznikiem (link function); jest to funkcja monoto-_, niczna i r´o˙zniczkowalna.

Tak wiec zale˙zno´s´c Y_, ize zmiennymi obja´sniajacymi nie jest bezpo´srednia, ale po´srednia,_, poprzez zastosowanie funkcji-linku g(.).

Powtarzajac: Teoria uog´olnionych modeli liniowych (GLM) wi_, a˙ze warto´sci oczeki-_, wane µ_i z obserwowanymi wektorami x_i poprzez funkcje g. Struktura zale˙zno´sciowa_, 10.12 mo˙ze by´c zapisana w innej postaci, uÃlatwiajacej obliczanie funkcji odwrotnej i_, pochodnej:

g(µ_i) = η_i, gdzie η_i = x^T_i β (10.13) Tak wiec mamy_,

Y₁, . . . , Y_N - obserwowane zmienne losowe, wzajemnie niezale˙zne µ₁, . . . , µ_N - warto´sci oczekiwane wyznaczone jako µ_i = E(Y_i) x₁, . . . , x_N - odpowiadajace im wektory zmiennych obja´sniaj_, acych_, η1, . . . , ηN - odpowiadajace im kombinacje liniowe η_, i = x^T_i β.

W ten spos´ob, za pomoca r´owna´n (10.12) i (10.13), – dla dowolnego osobnika i, dla_, kt´orego uzyskali´smy obserwacje Y_{, i} – zostaÃlo okre´slone, ˙ze :

warto´s´c oczekiwana zmiennej Y_i tego osobnika, czyli warto´s´c µ_i jest funkcja_, (kombinacja) liniow_, a warto´sci x_{, i1}, . . . , x_id stanowiacych zmienne obja´sniaj_, ace_, charakteryzujace tego osobnika._,

SkÃladowe wektora β = (β₁, . . . , β_d) okre´slajace kombinacj_, e liniow_, a x_, ^T_i β sa wsp´olne dla_, i = 1, . . . , N i stanowia parametery uog´olnionego modelu liniowego._,

Tym samym rozkÃlad obserwowanej warto´sci Y_i, i = 1, . . . , N , jest opisany funkcja_, gesto´sci postaci:_,

f_i(y) = f (y_i; x_i, β).

Je˙zeli funkcja f (.) nale˙zy do rodziny wykÃladniczej, to estymacja parametru β jest stosunkowo prosta. Szczeg´oÃlowe wzory – wraz z ich wyprowadzeniem – dla estymator´ow najwiekszej wiarygodno´sci mo˙zna znale´z´c w [3], str 39–41, lub [7], str 21–40._,

Ostateczny wynik jest taki, ˙ze estymatory te mo˙zna wyznaczy´c w spos´ob iteracyjny za pomoca iteracyjnej metody najmniejszych kwadrat´ow (Iterative Reweighted Least Squa-_, res). W kroku (k +1) tego algorytmu rozwiazuje si_, e ukÃlad r´owna´n liniowych, a rozwi_, azanie_, tego ukÃladu (oznaczane dalej jako β^(k+1)) dostarcza kolejnego przybli˙zenia dla wektora β dajacego maksimum funkcji wiarogodno´sci [3], str 39–41._,

W k + 1-tej iteracji rozwiazujemy nast_, epuj_, acy ukÃlad r´owna´n (symbole y_{, i} oznaczaja_, zaobserwowane realizacje zmiennych losowych Y_i): {irls }

X^TWXβ^(k+1) = X^TWz, (10.14)

gdzie

X_{N ×(d+1)} jest tablica danych poszerzon_, a o kolumn_, e jedynek,_,

β^(k+1) o wymiarach (d + 1) × 1 jest szukanym wektorem parametr´ow, otrzymamy go rozwiazuj_, ac ukÃlad r´owna´n (10.14),_,

z = z(k) = (z1, . . . , zN)^T jest wektorem pomocniczym, powstaÃlym z rozwiniecia 1-go_, rzedu (1st order) funkcji g(y) ≈ g(µ) + (y − µ)g_, ⁰(µ) = η + (y − µ)∂η/∂µ:

z_i = ˆη_i+ (y_i− ˆµ_i)∂ηi

∂µ_i =

d+1X

j=1

x_ijβ_j^(k)+ (y_i − µ_i)∂ηi

∂µ_i, z warto´sciami µ_i i ∂η_i/∂µ_i wyznaczonymi w punkcie β^(k),

W = W(k) jest macierza diagonaln_, a diag{w_{, ii}} ⁶ rozmiaru N × N, kt´orej elementy stanowia wariancj_, e rozwini_, ecia g(y) wyznaczon_, a przy zaÃlo˙zeniu. ˙ze η_{, i} i µ_i sa ustalone i_, znane (wyznaczone w poprzedniej, k-tej, iteracji):

w_ii⁻¹ = var(yi)(∂η_i

∂µ_i)²,

UkÃlad r´owna´n (10.14) przypomina ukÃlad r´owna´n normalnych dla modelu liniowego (np.

modelu regresji wielokrotnej), jednak – wobec faktu zale˙zno´sci W i z od aktualnej warto´sci β^(k) – musi by´c rozwiazywany iteracyjnie, startuj_, ac z jakiej´s warto´sci pocz_, atkowej β_, ⁽⁰⁾

Teoria GLM pozwala nam znale´z´c nie tylko estymatory ˆβ paramater´ow β, ale r´ownie˙z ich wariancje i kowariancje. Mamy, por. Dobson [3], str 63:

cov(ˆβ) = I⁻¹,

gdzie I jest warto´scia oczekiwan_, a hesjanu (macierz_, a informacji) o elementach_, I_hk= E[{∂²log L/∂β_h∂β_k}], h, k = 1, . . . , d + 1

wyznaczonego z logarytmu wiarygodno´sci L w punkcie ˆβ. Pokazuje sie, ˙ze I = X_, ^TWX.

6elementy wii oznaczaja tutaj wagi kolejnych wierszy tablicy danych X i nie maj_, a nic wsp´olnego z_, wagami sieci neuronowej

10.8 Metoda GLM – uog´olniony model liniowy 15

10.8.3 Zastosowanie GLM w klasyfikacji do 2 klas

Obserwowana zmienna losowa Yi ma rozkÃlad binarny (Bernoulliego), z prawdopodobie´n- stwem sukcesu π_i, gdy:

P r{Y_i} =

( 1, z p-stwem π_i, 0, z p-stwem 1 − π_i.

W zagadnieniach klasyfikacyjnych w przypadku klasyfikacji do dwu grup parametr π_i jest interpretowany jako oczekiwane prawdopodobie´nstwo sklasyfikowania osobnika nr i do klasy o etykiecie ’1’; natomiast wyra˙zenie 1 − πi oznacza zaklasyfikowanie osobnika nr i do klasy komplementarnej (oznaczonej umownie etykieta ’0’ lub ’2’)._,

Wska´znik ’i’ przy p-stwie π_i oznacza mo˙zliwo´s´c, ˙ze ka˙zdy osobnik (nr ’i’) m´ogÃl zosta´c wylosowany z rozkÃladu binarnego (Bernoulliego) o innym prawdopodobie´nstwie sukcesu.

Wiadomo, ˙ze w rozkÃladzie binarnym warto´s´c oczekiwana zmiennej losowej Y_ijest r´owna π_i. Oznaczmy og´olnie

µi = E(Yi) = πi.

Tak wiec ka˙zdy element ci_, agu {Y_{, i}} mo˙ze mie´c inna warto´s´c oczekiwan_, a. Warto´s´c ta mo˙ze_, zale˙ze´c od pewnych dodatkowych zmiennych, nazywanych zmiennymi obja´sniajacymi. Zmien-_, ne te sa dane, dla ka˙zdego elementu ci_, agu {Y_{, i}}, w postaci wektora x_i = (x_i1, . . . , x_id, 1)^T.

Niech η_i oznacza kombinacje liniow_, a rozwa˙zanych zmiennych obja´sniaj_, acych:_, η_i = x^T_i β, i = 1, . . . , N,

gdzie β oznacza wektor nieznanych wsp´oÃlczynnik´ow (w dalszym ciagu postaramy si_, e je_, wyestymowa´c na podstawie przyjetego modelu i zaobserwowanych warto´sci y_{, 1}, . . . , y_N).

Przyjmijmy, ˙ze warto´s´c oczekiwana Y_i, czyli parametr π_i zale˙zy tylko od kombinacji liniowej ηi = x^T_i β, a zale˙zno´s´c ta jest postaci {wz´or flogi}

E(Yi) = µi = πi = 1

1 + exp{− η_i} = 1

1 + exp{− x^T_i β} . (10.15) Chcac stosowa´c metody GLM, przyjmujemy jako funkcj_, e-Ãl_, acznik (link-function) funkcj_, e_, logit {wz´or logit}:

g(µ_i) = log( π_i

1 − π_i) = η_i = x^T_i β. (10.16) Poka˙zemy teraz, jak zdefiniowana funkcja logit zale˙zy od wektora zmiennych obja´snia- jacych x_{, i} i od wektora parametr´ow β. Mamy (M = [1 + exp(−η_i)], jest to mianownik we wzorze 10.15)

g(µi) = logit(µi) = log( πi

1 − π_i) = log( 1/M

exp(−η_i)/M) = ηi = x^T_i β.

Tak wiec wprowadzona logi speÃlnia warunki funkcji-Ãl_, acznika w uog´olnionym modelu linio-_, wym.

Poka˙zemy teraz, ˙ze funkcja odwrotn_, a do funkcji logit jest funkcja logistyczna._, Mamy bowiem

log( π_i

1 − π_i) = ηi ⇔ log(1 − π_i

π_i ) = − ηi ⇔ 1 − πi = πiexp{− ηi}.

Stad – wyznaczaj_, ac z ostatniej r´owno´sci powy˙zej – warto´s´c π_, i otrzymujemy πi = 1/[1 + exp{−η_i}], czyli wz´or (10.15).