)=E= @OIHOE=?O= IFAJ @ = 5J=JOIJO?A AJ@O ==EO @=O?D

Analiza dyskryminacyjna

Konspekt do zaj¦¢: Statystyczne metody analizy danych

Agnieszka Nowak-Brzezi«ska 8 stycznia 2010

1 Wprowadzenie

Celem zaj¦¢ ma by¢ omówienie podstawowych zagadnie« zwi¡zanych z analiz¡

dyskryminacyjn¡. Liniowa Analiza Dyskryminacyjna i zwi¡zana z ni¡ Fishe- rowska liniowa analiza dyskryminacyjna s¡ u»ywane w uczeniu maszynowym do znalezienia liniowej kombinacji cech, które najlepiej rozró»niaj¡ dwie lub wi¦cej klas obiektów lub zdarze«. Wynikowe kombinacje s¡ u»ywane jako klasykator liniowy lub, cz¦±ciej, sªu»¡ redukcji wymiarów do pó¹niejszej klasykacji sta- tystycznej. LDA jest u»ywana do wielu zastosowa« zwi¡zanych z klasykacj¡.

Jednym z nich jest rozpoznawanie twarzy. Obraz twarzy skªadaj¡cy si¦ z bardzo du»ej ilo±ci pikseli jest redukowany do mniejszego zbioru linowych kombinacji, które mog¡ nast¦pnie by¢ wykorzystane do klasykacji.

Tradycyjnie ju», podstawowym ¹ródªem informacji w tym zakresie powinien by¢ podr¦cznik: [ J. Koronacki, J. ‚wik: Statystyczne systemy ucz¡ce si¦, wy- danie drugie, Exit, Warsaw, 2008. (rozdziaª 1)].

2 Analiza dyskryminacyjna

Gdyby±my zerkn¦li do dost¦pnych ¹ródeª po denicj¦ poj¦cia dyskryminacji otrzymamy nast¦puj¡c¡ odpowied¹:

Dyskryminacja (ªac. discriminatio - 'rozró»nianie'; czyt.: diskri- minacjo) - cz¦sto wyst¦puj¡ca forma wykluczenia spoªecznego, obja- wiaj¡ca si¦ poprzez traktowanie danej osoby mniej przychylnie, ni»

innej w porównywalnej sytuacji ze wzgl¦du na jak¡± cech¦ (np. pªe¢, to»samo±¢ seksualna, wiek, niepeªnosprawno±¢, religia lub przekona- nia czy pochodzenie etniczne lub rasowe).

Odnosz¡c t¦ denicj¦ do obszaru metod statystycznych, znajdziemy metod¦

nazywan¡ analiz¡ dyskryminacyjn¡. Metod¡ u»ywan¡ do znalezienia liniowej kombinacji cech, które najlepiej rozró»niaj¡ dwie lub wi¦cej klas obiektów lub zdarze« jest liniowa analiza dyskryminacyjna (ang. Linear discriminant analy- sis, LDA). Okre±la si¦ j¡ mianem metody klasykacji pod nadzorem.

Analiza dyskryminacyjna d¡»y do tego, aby ka»dej obserwacji x ∈ X przy- pisa¢ klas¦, do której ta obserwacja nale»y. Reguªa dyskryminacyjna (klasyka- cyjna) dzieli wi¦c zbiór X na g rozª¡cznych podzbiorów (klas) i daje si¦ zapisa¢

jako odwzorowanie: d(x) : X → G, gdzie G jest sko«czonym g-elelemntowym zbiorem etykiet klas, do których nale»¡ obserwacje, przy czym g ≥ 2.

2.1 Nierozª¡czno±¢ klas

Pami¦tajmy, »e typowe procesy klasykacyjne obarczone s¡ zwykle niepewno-

±ci¡. Zawsze mo»emy w próbie ucz¡cej znale¹¢ obserwacje nietypowe, które wprowadz¡ szum informacyjny i zaburz¡ nasz¡ klasykacj¦. Nie wolno nam zakªada¢ rozdzielno±ci klas cho¢ oczywi±cie tego by±my sobie najch¦tniej »y- czyli buduj¡c reguªy klasykacyjne. Fakt konieczno±ci wpuszczenia niepewno±ci do takich modeli mo»e by¢ wspomagany probabilistyk¡, i wymaga traktowania obserwacji w próbie ucz¡cej jako wektorów losowych.

3 Fisherowska analiza dyskryminacyjna

Analiza dyskryminacyjna w uj¦ciu Fishera, jest okre±lana jako historycznie pierwsze podej±cie do klasykacji pod nadzorem. Zakªada, »e wektory obserwa- cji s¡ wektorami w przestrzeni p - wymiarowej X ⊂ R^pi prowadzi do otrzymania reguªy dyskryminacyjnej opartej na funkcji liniowej. Reguªa ta dla g = 2 stara si¦ znale¹¢ kierunek a w X, który najlepiej rozdziela obydwie podpróby ucz¡ce, przy tym konstrukuj¡c miar¦ odlegªo±ci mi¦dzy klasami, uwzgl¦dni¢ zmienno±¢

wewn¡trzgrupow¡. Jasne jest, »e rozproszenie wewn¡trzgrupowe powinno by¢

scharakteryzowane odpowiednimi - opartymi na stosownych podpróbach - prób- kowymi macierzami kowariancji. Konstrukcja metody LDA wg Fishera wymaga zagregowania informacji o klasach, do których mamy klasykowa¢ nowe obser- wacje, za pomoc¡ wska¹ników poªo»enia i rozproszenia g podprób próby ucz¡cej.

Je±li mamy tylko dwie klasy, to prób¦ (x1, y1), (x2, y2), . . . , (xn, yn)mo»emy po- dzieli¢ na podprób¦ obserwacji z klasy pierwszej i podprób¦ obserwacji z klasy drugiej:

x11, x12, . . . , x1n1 z klasy 1, x21, x22, . . . , x2n₂ z klasy 2,

gdzie n = n1+ n2. Wtedy ±rednie klas mo»emy zapisa¢ jako:

xk= 1 nk

n_k

X

i=1

xki dla k = 1, 2.

Je±li jako xkokre±limy ±redni¡ grupow¡ i zaªo»ymy, »e wariancja wektora X b¦dzie taka sama we wszystkich populacjach (tutaj w dwóch, bo k = 1, 2), to wspólna macierz kowariancji zostanie wyznaczona poprzez obliczenie macierzy kowariancji wewn¡trzgrupowej W :

W = 1

n − 2 X2 k=1

(nk− 1)Sk= 1 n − 2

X2 k=1

n_k

X

l=1

(xkl− xk)(xkl− xk)⁰

gdzie Sk oznacza macierz kowariancji próby w k-tej populacji, x⁰ -transpozycj¦

wektora x, i powiemy tak»e, »e n = n1+ n2. Fisher zadanie analizy dyskrymi- nacyjnej zdeniowaª nast¦puj¡co:

znale¹¢ taki kierunek a, który maksymalizuje odlegªo±¢ mi¦dzy zrzu- towanymi ±rednimi obu prób przy uwzgl¦dnieniu wariancji rzutu

(czyli standaryzowanej odlegªo±ci zrzutu ±rednich)

argmaxa(q a⁰x2− a⁰x1

a⁰W a(_n¹

1 +_n¹

2))²= argmaxa(a⁰x2− a⁰x1)² a⁰W a Rozwi¡zaniem jest:

a = W⁻¹(x2− x1)

ale samo wyznaczenie a nie tworzy jeszcze reguªy dyskryminacyjnej. By j¡

zbudowa¢:

ˆ rzutujemy obserwacje x na kierunek a,

ˆ i nast¦pnie klasykujemy je do pierwszej lub drugiej klasy, w zale»no±ci od tego czy rzut ten byª bli»szy rzutowi ±rodka próby pierwszej czy drugiej:

d = 2 ⇔ a⁰(x − (x1+ x2)/2) > 0.

Idea analizy dyskryminacyjnej zostaªa bardzo dobrze oddana na ilustracji w podr¦czniku Prof. Koronackiego pt. Statystyczne systemy ucz¡ce si¦ w rozdziale 1 na stronie 30. Zrzut tej ilustracji przedstawia rysunek 1.

Rysunek 1: Wykres obserwacji dla klas 1 (trójk¡tów) i 2 (kóªek)

Widzimy tutaj, »e obserwacje z dwóch klas (kóªek i trójk¡tów) s¡ prezentowane w przestrzeni dwuwymiarowej jako punkty. Punkt x1reprezentuje ±redni¡ grupy 1, za± x2 odpowiednio ±redni¡ grupy drugiej. Na rysunku przedstawiono rzuty obserwacji i ±rednich z prób dla kazdej klasy na prost¡ opisan¡ jako ea. Do tej prostej poprowadzono prost¡ - nazwan¡ na wykresie prost¡ dyskryminacyjn¡.

3.1 Reguªa dyskryminacyjna

Reguªa dyskryminacyjna wyznacza w zbiorze X hiperpowierzchni¦ dyskrymina- cyjn¡ dziel¡c¡ ten zbiór na dwa podzbiory: tych punktów x, które nale»¡ do klasy 1 i tych, które nale»¡ do klasy 2. Co wa»ne:

reguªa taka wyznacza hiperpowierzchni¦, która jest prosto- padªa do kierunku ea i przechodzi przez ±rodek odcinka ª¡- cz¡cego punkty ea⁰x1 oraz ea⁰x2.

Je±li spojrzymy na wszystkie punkty le»¡ce po tej samej stronie co x1 to oka»e si¦, »e przecie» i ich rzuty ortogonalne na kierunek ea musz¡ le»e¢ po tej samej stronie, a wi¦c speªniaj¡ warunek: |ea⁰x − ea⁰x₁| < |ea⁰x − ea⁰x₂| i b¦d¡

klasykowane do klasy 1.

Reguªa dyskryminacyjna Fishera Reguªa ta mówi¢ b¦dzie, »e

przypiszemy obiekt x ∈ X do klasy 2 je±li (x2− x1)⁰W⁻¹[x −¹₂(x1+ x2)] > 0oraz do klasy 1 w przypadku, gdy prawdziwa jest nierówno±¢

odwrotna. Kontrowersyjny przypadek, gdy punktu x le»y dokªad- nie na hiperpªaszczy¹nie dyskryminacyjnej rozwi¡zuje si¦ poprzez zrzucenie decyzji na eksperta.

4 Problem wielu klas

Koncepcj¦ Fishera z liczb¡ klas g równ¡ 2 mo»a uogólni¢ na przypadek wi¦kszej liczby klas (czyli gdy g > 2). Takie uogólnienie jako pierwsi zaproponowali C.R.Rao i J.G.Bryan.

Wtedy równie» b¦dziemy szuka¢ takiego kierunku a, który pozwoli maksy- malizowa¢ wyra»enie _a^a0⁰W a^Ba, gdzie:

B = 1 g − 1

Xg k=1

nk(xk− x)(xk− x)⁰

jest macierz¡ kowariancji mi¦dzygrupowej za±

W = 1

n − g Xg k=1

(nk− 1)Sk

jest macierz¡ kowariancji wewn¡trzgrupowej.

4.1 Reguªa dyskryminacyjna dla wielu klas

Maj¡c szukany wektor ea obserwacj¦ x przypiszemy dla j-tej klasy, gdy |ea⁰x − e

a⁰xj| < |ea⁰x − ea⁰xk| dla wszystkich j 6= k. Najcz¦±ciej nie jest mo»liwe okre-

±lenie jednego kierunku rozdzielaj¡cego klasy obserwacji. Wtedy szuka si¦ np.

dwóch kierunków kanonicznych. Z tego wzgl¦du - lepiej rozwi¡za¢ ( g

2 )zada«

dyskryminacyjnych mi¦dzy dwiema klasami, co po prostu oznacza rozwi¡zanie

zadania dyskryminacji dla wszystkich mo»liwych zestawów par klas spo±ród g klas.

Uogólnienie reguªy Fishera dla dwóch klas na reguª¦:

zaklasykuj wektor x do populacji j, je±li |a⁰x − a⁰xj| ≤ |a⁰x − a⁰xk| dla k 6= j, gdzie a jest pierwszym wektorem kanonicznym

nie zawsze sprawdza si¦ dla wi¦cej ni» dwóch populacji.

Dla wszystkich zestawów par klas budujemy wspóln¡ macierz kowariancji wewn¡trzgrupowej. Podziaªu przestrzeni R^p dokonuj¡ hiperpªaszczyzny dys- kryminacyjne otrzymane z reguªy klasykacyjnej dij skonstruowanej dla dowol- nej pary i, j populacji. Co wa»ne, te hiperpªaszczyzny nie musz¡ si¦ przecina¢

w jednym punkcie. Je±li rozwi¡zujemy zadanie dyskryminacji dla wszystkich mo»liwych zestawów par klas spo±ród g klas, wówczas dla tych wszystkich klas przyjmujemy wspóln¡ macierz kowariancji wewn¡trzgrupowej. Gdyby±my ma- cierze obliczali oddzielnie dla pary klas 1 i 2, pary klas 1 i 3 oraz pary klas 2 i 3, otrzymane na tej podstawie proste dyskryminacyjne mogªyby si¦ nie przeci¡¢ w jednym punkcie - ale tworzyªyby trójk¡t, w którym nie byªoby jednoznacznego rozwi¡zania zadania dyskryminacji.

4.2 Okre±lenia stosowane w analizie dyskryminacyjnej

Zmienn¡ ea⁰x nazywa¢ b¦dziemy pierwsz¡ zmienn¡ kanoniczn¡ odpowiadaj¡c¡

wektorowi x, a wektor ea pierwszym wektorem kanonicznym. ea jest kierunkiem kanonicznym uzyskanym przez maksymalizacj¦ ilorazu _a^a00_{W a}^Ba i jest to takie kie- runek, który najlepiej rozdziela klasy.

4.2.1 Wektory kanoniczne i zmienne kanoniczne

Maj¡c dan¡ macierz W⁻¹B, gdzie W jest macierz¡ kowariancji wewn¡trzgru- powej, za± B kowariancji mi¦dzygrupowej, powiemy, »e kolejne wektory wªasne odpowiadaj¡ce uporz¡dkowanym od najwi¦kszej do najmniejszej warto±ciom wªasnym tej macierzy nosz¡ nazw¦ wektorów kanonicznych. Wektor a zde- niowany w metodzie Fishera jest pierwszym wektorem kanonicznym dla g = 2.

Wektorów kanonicznych jest maksymalnie g − 1. Je±li ai jest i-tym wektorem kanonicznym, to a⁰ixnazywana jest i-t¡ zmienn¡ kanoniczn¡. W problemach wielowymiarowych u»ywa si¦ cz¦sto dwoch lub trzech pierwszych zmiennych ka- nonicznych do wizualizacji wektora cech. Zmienne kanoniczne daj¡ce wyra¹n¡

separacj¦ grup na wykresie s¡ z reguªy istotne w klasykacji.

4.2.2 Liczba kierunków kanonicznych

Przy wi¦cej ni» dwóch klasach, gdzie celem jest znalezienie nie jednego ale wielu kierunków kanonicznych powstaje problem optymalnego oszacowaniach ich liczby. Wa»nym jest tak»e rozwi¡zanie zadania wielokrotnej maksymalizacji

a⁰Ba

a⁰W a, które nierzadko przynosi bardzo ciekawe uproszczenie oryginalnego pro- blemu przez redukcj¦ jego wymiaru i rozwi¡zanie zadania dyskryminacji (kla- sykacji) w otrzymanej przestrzeni o zmniejszonym wymiarze, a mianowicie w t-wymiarowej przestrzeni euklidesowej t zmiennych kanonicznych. Zwykle wy- miar p przestrzeni obserwacji jest wi¦kszy, lub znacznie wi¦kszy, od liczby klas

pomniejszonej o 1, czyli g − 1. Je±li zast¡pimy oryginalny problem wielowy- miarowy problemem o znacznie mniejszym wymiarze równym g − 1 lub nawet mniejszym, za to w przestrzeni generowanej przez kierunki najlepiej rozdziela- j¡ce klasy, uzyskamy w prosty sposób zadowalaj¡ce wyniki klasykacji. Wtedy bowiem zadanie klasykacji mo»emy wykona¢ stosuj¡c dowoln¡ z metod po- wszechnie znanych i stosowanych w problemach klasykacji.

5 Analiza dyskryminacyjna za pomoc¡ R

Do wykonania ¢wicze« wykorzystamy przykªadowe dane udost¦pnione przez au- torów podr¦cznika:

J. ‚wik, J. Mielniczuk: Statystyczne systemy ucz¡ce si¦ - ¢wiczenia w oparciu o pakiet R, Ocyna Wydawnicza PW, Warszawa, 2009.

Zródªo pobrania danych to:

http://www.ipipan.waw.pl/staff/j.mielniczuk/SSUS-Programy-Dane.zip

5.1 Dwie klasy i dwie zmienne

Jako pierwszy do analizy we¹miemy prosty zbiór earthquake dotycz¡cy klasy-

kacji wybuchów nuklearnych i trz¦sie« ziemi. Zmienn¡ klasykuj¡c¡ jest tutaj zmienna popn, która mo»e przyjmowa¢ jedn¡ z dwóch mo»liwych warto±ci: equ- ake lub explosn. Do pierwszej klasy odpowiednio zaliczymy 20 obserwacji, do drugiej - 9. Wczytanie zbioru danych mo»liwe jest dzi¦ki komendzie:

> earthquake<-read.table("C:\\earthquake.txt",header=TRUE)

Je±li chcemy by obserwacje z ró»nych klas byªy oznaczone innymi kolorami, umo»liwi to komenda:

> kolor <-c(rep("black",20),rep("red",9))

Teraz je±li chcemy tak»e by byªy takie obserwacje z ró»nych klas etykietowane innymi symbolami, wystarczy u»y¢ komendy:

> symbole <-c(rep("Q",20),rep("X",9)) Teraz stosuj¡c instrukcj¦

> plot(earthquake$body, earthquake$surface, type="n",xlab="body",ylab="surface")

> text(earthquake$body,earthquake$surface,symbole,col=kolor)

w efekcie uzyskamy wykres rozproszenia obiektów zbioru earthquake dla zmiennych body oraz surface (rysunek 2).

5.2 Wyznaczenie równania prostej rozdzielaj¡cej dwie klasy z denicji

Kolejno musimy wykona¢ kroki:

ˆ okre±li¢ liczno±ci próby i klas:

4.5 5.0 5.5 6.0 6.5

4.04.55.05.56.0

body

surface

Q

Q

Q

Q

Q

Q Q

Q

Q

Q Q

Q

Q Q Q

QQ

Q Q

Q

XX X X

X

X X

X X

Rysunek 2: Wykres rozproszenia obiektów zbioru earthquake

> n<-29

> n1<-20

> n2<-9

ˆ wyznaczy¢ macierze kowariankcji - które zakªadamy, »e s¡ równe:

> kowar1 <- cov(earthquake[1:20,2:3])

> kowar2 <- cov(earthquake[21:29,2:3])

ˆ wyznaczy¢ macierz kowariancji wewn¡trzgrupowej:

> kowar <- ((n1-1)*kowar1+ (n2-1)*kowar2)/(n-2)

ˆ potrzebujemy te» macierzy odwrotnej do macierzy kowar. U»yjemy do tego funkcji ginv z pakietu MASS:

> library(MASS)

> invkowar<-ginv(kowar)

ˆ musimy jeszcze okre±li¢ wektory ±rednich grupowych osobno dla obserwacji z ró»nych grup (wiadomo, »e obserwacje o indeksach 1..20 nale»¡ do jednej klasy, a te o indeksach 21..29 do drugiej klasy):

> sr1<-mean(earthquake[1:20,2:3])

> sr2<-mean(earthquake[21:29,2:3])

ˆ na koniec nale»y okre±li¢ wspóªczynniki a1 i a2 oraz wyraz wolny prostej x2 = −(a1/a2) ∗ x1 − stala/a2rozdzielaj¡cej klasy:

> wspa <-invkowar %*% (sr2-sr1)

> stala<-0.5&(t(sr1+sr2) %*% wspa)

ˆ wyrysowanie takiej prostej wykonamy komend¡:

> abline(-stala/wspa[2],-wspa[1]/wspa[2])

5.3 Analiza dyskryminacyjna z u»yciem funkcji lda

W pakiecie R do analizy dyskryminacyjnej mo»na wykorzysta¢ pakiet lda {MASS}.

Jej u»ycie dla zbioru earthquake wygl¡da¢ b¦dzie nast¦puj¡co:

> equake.lda <-lda(popn~.,data=earthquake)

> print(equake.lda) Call:

lda(popn ~ ., data = earthquake) Prior probabilities of groups:

equake explosn 0.6896552 0.3103448 Group means:

body surface equake 5.249000 4.740000 explosn 5.978889 4.244444

Coefficients of linear discriminants:

LD1 body 3.918916 surface -1.865219

Je±li wi¦c estymatory prawdopodobie«stw przynale»no±ci klasowej wynosz¡

odpowiednio:

equake explosn 0.6896552 0.3103448

to równanie prostej b¦dzie miaªo posta¢ komendy:

> abline((-stala +log(0.3103448/0.6896552))/wspa[2],-wspa[1]/wspa[2],lty=2)

5.4 Reklasykacja

Tabela reklasykacji przy u»yciu LDA 29 obiektów o znanej przynale»no±ci do klasy przedstawia si¦ nast¦puj¡co:

> equake.pred<-predict(equake.lda,newdata=earthquake)

> print(table(earthquake$popn,equake.pred$class)) equake explosn

equake 19 1

explosn 0 9

Widzimy tu, »e jeden z obiektów, który nale»aª do klasy equake zostaª bª¦d- nie wrzucony do klasy explosn. Natomiast w przypadku drugiej klasy bª¦dów klasykacji nie zauwa»ono. Zªa klasykacja jest zobrazowana jako takie punkty na wykresie, które s¡ np. blisko linii przerywanej. Mo»emy wyznaczy¢, które punkty b¦d¡ ¹le sklasykowane poprzez obliczenie prawdopodobie«stw aposte- riori zaklasykowania do ka»dej z klas. Wywoªanie komend:

> print(equake.pred$post) equake explosn 1 3.307965e-01 6.692035e-01 2 9.632382e-01 3.676183e-02

3 9.621776e-01 3.782240e-02 4 9.895360e-01 1.046398e-02 5 9.998841e-01 1.158851e-04 ....

27 4.775093e-02 9.522491e-01 28 5.307004e-04 9.994693e-01 29 6.599999e-04 9.993400e-01

Widzimy tutaj, »e dla obserwacji pierwszej prawdopodobie«stwo zaklasy- kowania do klasy 1 wynosi 3.307 za± do klasy 2: 6.692 co oznacza, »e po prostu prawdopodobienstwo zaklasykowania tego obiektu do klasy 2 jest wi¦ksze ni»

do klasy 1 - st¡d bª¡d klasykacji.

5.5 Klasykacja z u»yciem zmiennych kanonicznych

Dla dwóch klas (g = 2) otrzymamy jeden wektor kanoniczny a1i jedn¡ zmienn¡

kanoniczn¡ a⁰1x. Wyznaczenie warto±ci wektora kanonicznego umo»liwia ko- menda R:

> print(equake.lda$scaling) czego efektem b¦d¡ warto±ci:

body 3.918916LD1 surface -1.865219

Z kolei wyznaczenie warto±ci pierwszej zmiennej kanonicznej a⁰1xdla elementów próby b¦dzie mo»liwe dzi¦ki komendzie:

> kanon<-t(equake.lda$scaling)%*%t(as.matrix(earthquake[,2:3]))

> print(kanon)

której efektem b¦d¡ dane:

[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]

LD1 14.01875 12.96967 12.97748 12.63056 11.43799 12.36923 10.74818 10.40789

[,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]

LD1 10.72615 11.27759 12.42707 11.90481 12.6895 13.01452 11.28701 10.80771

[,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]

LD1 11.19206 10.26016 10.38040 11.05831 15.59385 15.20762 14.77277 14.96852

[,25] [,26] [,27] [,28] [,29]

LD1 15.72968 17.13731 14.62335 15.82502 15.76738

Progiem rozdzielaj¡cym klasy w podejsciu Fishera jest warto±¢: a⁰(x1+ x2)/2. Je±li ±rednie grupowe wyznaczymy komend¡ > print(equake.lda\$means)i b¦d¡ one nast¦puj¡ce:

body surface equake 5.249000 4.740000 explosn 5.978889 4.244444

to próg rozdzielaj¡cy klasy popn i equake ma warto±¢ wyznaczan¡ jako:

> print(3.918916*(0.5*(5.249000 + 5.978889))-1.865219*(0.5*(4.74+4.244444)))

która wynosi: [1] 13.6216. Jak wida¢ byªa ona wyznaczona z równiania:

LD1body∗((equakebody+explosnbody)/2)+LD1surf ace∗((equakesurf ace+explosnsurf ace)/2) Wtedy reguªa klasykacyjna ma posta¢ nast¦puj¡c¡:

ˆ je±li warto±¢ zmiennej kanonicznej dla jakiej± obserwacji jest mniejsza od warto±ci 13.6216 to t¦ obserwacj¦ zaliczymy do klasy pierwszej (equake),

ˆ w przeciwnym przypadku - do klasy drugiej (explosn).

Patrz¡c teraz na warto±ci zmiennej kanonicznej kanon mo»emy sprawdzi¢, czy s¡ takie obserwacje, które powinny na podstawie warto±ci funkcji kanonicznej nale»e¢ do klasy pierwszej a zostaªy niepoprawnie przydzielone do klasy drugiej.

Tak jest w przypadku obserwacji pierwszej. Pami¦tamy bowiem, »e pierwsze 20 obserwacji to klasa equake. Je±li jednak b¦dziemy klasykowa¢ obserwacje wg warto±ci kanonicznych, to warto±¢ tej obserwacji równa 14.01875 sugeruje, »e obserwacja ta nale»y do klasy explosn skoro ma warto±¢ wi¦ksz¡ ni» progowa warto±¢ 13.6216. Wykres warto±ci pierwszej zmiennej kanonicznej dla elementów próby przedstawia rysunek 3. Do jego wygenerowania niezb¦dne jest wygene- rowanie nast¦puj¡cego kodu:

> equake.pred<-predict(equake.lda,newdata=earthquake)

> plot(equake.pred$x,type="n",xlab="index",ylab="pierwsza zmienna kanoniczna")

> symbole <-c(rep("Q",20),rep("X",9))

> text(equake.pred$x,symbole,cex=0.8)

> abline(0,0)

0 5 10 15 20 25 30

−2−101234

index

pierwsza zmienna kanoniczna

Q

Q Q Q

Q Q

Q Q

Q Q

Q Q

Q Q

Q Q

Q

QQ Q

X X

XX X

X

X X X

Rysunek 3: Wykres warto±ci pierwszej zmiennej kanonicznej

Wida¢ na tym wykresie, »e faktycznie je±li prosta rozdzialaj¡ca obie klasy b¦dzie w tym miejscu w którym jest teraz, wówczas obserwacja numer 1 zostanie ¹le zaklasykowana.

6 3 klasy, 2 zmienne obja±niaj¡ce

Dobrze znanym i cz¦sto analizowanym zbiorem danych jest zbiór irysów, dziel¡cy 150 obserwacji do 3 klas: setosa, virginica i versicolor.

> irysy<-read.table("d:\\iris1.txt",header=TRUE)

> irysy.lda=lda(Species~Petal.Length+Petal.Width,data=iris)

> irysy.pred=predict(irysy.lda)

> print(tabl<-table(irysy$Species,irysy.pred$class))

Efektem jest wy±wietlenie poprawnie sklasykowanych elementów próby:

setosa versicolor virginica

setosa 50 0 0

versicolor 0 48 2

virginica 0 4 46

Je±li chcemy pozna¢ procent poprawnej klasykacji dla próby treningowej musimy jedynie u»y¢ komend:

> procent<-(100*sum(diag(tabl))/sum(tabl))

> print(procent)

W efekcie otrzymamy warto±¢:

[1] "Procent poprawnej klasyfikacji dla proby treningowej"

[1] 96

Wykres rozproszenia z zaznaczeniem numerów klas wyznaczymy za pomoc¡ ko- mendy:

> plot(irysy[,c(2,5)],type="n",xlab="Petal.Length",ylab="Petal.Width")

> text(irysy$Petal.Length,irysy$Petal.Width,as.character(irysy$Species),cex=0.8)

6.1 Dyskryminacja metod¡ QDA

Chc¡c u»y¢ kwadratowej funkcji dyskryminacyjnej dla tego samego zbioru obserwacji z trzema klasami musimy po prostu wywoªa¢ zamiast funkcji lda funkcj¦ qda:

> irysy.qda=qda(Species~Petal.Length+Petal.Width,data=iris)

> irysy.predQ = predict(irysy.qda)

> print(tabl<-table(irysy$Species,irysy.predQ$class)) W efekcie tabela reklasykacji przedstawia si¦ nast¦puj¡co:

setosa versicolor virginica

setosa 50 0 0

versicolor 0 49 1

virginica 0 2 48

Analoglicznie jak dla lda procent poprawnej klasykacji wyznaczymy poprzez wywo- ªanie komend:

> procent<-100*(sum(diag(tabl))/sum(tabl))

> print(procent)

Dla tego zbioru wynosi on [1] 98.

6.2 LDA dla zbioru iris dla wszystkich 4 zmiennych obja-

±niaj¡cych

Wczytanie caªego zbioru do analizy oraz wywoªanie funkcji lda wygl¡da nast¦puj¡co:

> irysy<-read.table("d:\\iris1.txt",header=TRUE)

> irysy.lda <-lda(Species~.,data=iris)

> print(irysy.lda)

W efekcie otrzymyjemy raport:

Call:

lda(Species ~ ., data = iris) Prior probabilities of groups:

setosa versicolor virginica 0.3333333 0.3333333 0.3333333 Group means:

Sepal.Length Sepal.Width Petal.Length Petal.Width

setosa 5.006 3.428 1.462 0.246

versicolor 5.936 2.770 4.260 1.326

virginica 6.588 2.974 5.552 2.026

Coefficients of linear discriminants:

LD1 LD2

Sepal.Length 0.8293776 0.02410215 Sepal.Width 1.5344731 2.16452123 Petal.Length -2.2012117 -0.93192121 Petal.Width -2.8104603 2.83918785 Proportion of trace:

LD1 LD2 0.9912 0.0088

>

Je±li chcemy wyrysowa¢ po prostu obserwacje z trzech ró»nych klas innymi kolorem i symbolem na wykresie to mo»liwe to b¦dzie przez wykonanie nast¦- puj¡cych instrukcji:

> kolor <-c(rep("black",50),rep("red",50),rep("blue",50))

> symbole <-c(rep("Q",50),rep("X",50),rep("Y",50))

> plot(irysy$Sepal.Width, irysy$Sepal.Length, type="n",xlab="",ylab="")

> text(irysy$Sepal.Width,irysy$Sepal.Length,symbole,col=kolor) Efektem b¦dzie wykres 4.

6.3 Inne fragmenty analizy z u»yciem R dla zbioru iris

Kod programu w ±rodowisku R, który zliczy liczebno±¢ klas ze zbioru iris, gdzie obserwacje z ró»nych klas maj¡ inne symbole (s, c lub v), i tylko dla losowych 75 obserwacji, wygl¡da nast¦puj¡co:

2.0 2.5 3.0 3.5 4.0

4.55.05.56.06.57.07.58.0

Q Q

QQ Q

Q

Q Q

Q Q

Q

Q Q

Q

Q Q

Q Q

Q

Q Q

Q

Q Q

Q

Q Q

Q Q

QQ Q

Q Q

QQ Q

Q

Q QQ

Q Q

Q Q

Q

Q

Q Q Q X

X X

X X

X X

X X

X X

X XX

X X

X X X

X

X X X

X X

X X X

X X X X

X X

X X X

X

X X X X X

X

X XX X

X X

Y

Y Y

Y Y Y

Y Y

Y

Y

Y Y Y

Y Y

Y Y

Y Y

Y

Y

Y Y

Y

Y Y

Y Y Y

Y Y

Y

YY Y

Y

Y Y Y

Y Y Y

Y

YY Y

Y Y

Y Y

Rysunek 4: Wykres rozrzutu obserwacji ze zbioru Iris

> data(iris3)

> Iris <- data.frame(rbind(iris3[,,1], iris3[,,2], iris3[,,3]), Sp = rep(c("s","c","v"), rep(50,3)))

> train <- sample(1:150, 75)

> table(Iris$Sp[train]) Efektem b¦dzie c s v

28 20 27

co oznacza jak ªatwo si¦ domy±le¢, »e w zbiorze ucz¡cym wyst¡piªo 28 obser- wacje z klasy oznaczonej c, 20 obserwacji z klasy s oraz 27 z klasy v. Je±li rozszerzymy t¦ analiz¦ na caªy zbiór obserwacji, który jak wiadomo w oryginale liczy 150 obserwacji efekt b¦dzie nast¦puj¡cy:

> train <- sample(1:150, 150)

> table(Iris$Sp[train]) c s v

50 50 50

Je±li teraz chcemy znale¹¢ dla takiego zbioru reguª¦ klasykacyjn¡, u»yjemy funkcji lda nast¦puj¡co:

> train <- sample(1:150, 75)

> table(Iris$Sp[train]) c s v

29 22 24

> z <- lda(Sp ~ ., Iris, prior = c(1,1,1)/3, subset = train)

> predict(z, Iris[-train, ])$class

[1] s s s s s s s s s s s s s s s s s s s s s s s s s s s s c c c c c v c c c c [39] c c c c c c c c c c c v v v v v v v v v v v v v v v v v v v v v v v v v v Levels: c s v

Zastosowali±my tu tak»e równocze±nie funkcj¦ sªu»¡c¡ predykcji, która dla ka»dej z 75 obserwacji podaªa przewidywan¡ klas¦ przynale»no±ci. Je±li np nie chcemy u»ywa¢ pewnej zmiennej w predykcji np. zmiennej Petal.W mo»emy to zrobi¢ nast¦puj¡co:

> (z1 <- update(z, . ~ . - Petal.W.)) Wówczas na ekranie zobaczymy raport:

Call:

lda(Sp ~ Sepal.L. + Sepal.W. + Petal.L., data = Iris, prior = c(1, 1, 1)/3, subset = train)

Prior probabilities of groups:

c s v

0.3333333 0.3333333 0.3333333 Group means:

Sepal.L. Sepal.W. Petal.L.

c 5.872414 2.713793 4.231034 s 5.009091 3.450000 1.477273 v 6.616667 2.983333 5.579167

Coefficients of linear discriminants:

LD1 LD2

Sepal.L. 1.293012 -0.4744279 Sepal.W. 1.342522 2.9403952 Petal.L. -3.250248 0.6057279 Proportion of trace:

LD1 LD2 0.9891 0.0109

>

W raporcie tym, widzimy prawdopodobie«stwa przynale»no±ci do poszcze- gólnych klas, ±rednie grupowe dla ka»dej klasy, przy uwzgl¦dnieniu oczywi±cie tylko tych trzech zmiennych obja±niaj¡cych:

Sepal.L., Sepal.W. oraz Petal.L..

komenda do wyrysowania takiego modelu analizy dyskryminacyjnej przechowa- nej w z wygl¡da nast¦puj¡co:

> plot(z, dimen=1)

a jej efektem b¦dzie rysunek 5.

Efektem wywoªania nast¦puj¡cego fragmentu programu w R:

−10 −5 0 5

0.00.20.4

group c

−10 −5 0 5

0.00.20.4

group s

−10 −5 0 5

0.00.20.4

group v

Rysunek 5: Analiza dyskryminacyjna dla zbioru Iris

> iris$klasa<-factor(rep(c("s","c","v"),rep(50,3)))

> attach(iris)

> lda.iris<-lda(klasa~Petal.Width+Petal.Length,data=iris)

oraz komendy > print(lda.iris) b¦dzie wy±wietlenie nast¦puj¡cych wy- ników:

Call:

lda(klasa ~ Petal.Width + Petal.Length, data = iris) Prior probabilities of groups:

c s v

0.3333333 0.3333333 0.3333333 Group means:

Petal.Width Petal.Length

c 1.326 4.260

s 0.246 1.462

v 2.026 5.552

Coefficients of linear discriminants:

LD1 LD2

Petal.Width -2.402394 5.042599 Petal.Length -1.544371 -2.161222 Proportion of trace:

LD1 LD2 0.9947 0.0053

Taki raport najpierw przedstawia nam prawdopodobie«stwa dla klas obli- czone na podstawie zbioru ucz¡cego. Poniewa» wiadomo, »e w zbiorze iris 150 obserwacji jest równo podzielone na 3 klasy, w ka»dej po 50 obiektów, tak

te» równo rozkªadaj¡ si¦ tu prawdopodobie«stwa. Cz¦±¢ raportu mówi¡ca o wspóªrz¦dnych ±rodków ci¦»ko±ci to blok okre±lony jako:

Group means:

Petal.Width Petal.Length

c 1.326 4.260

s 0.246 1.462

v 2.026 5.552

Liniowe funkcje dyskryminacyjne odczytamy z kolejnego fragmentu raportu:

Coefficients of linear discriminants:

LD1 LD2

Petal.Width -2.402394 5.042599 Petal.Length -1.544371 -2.161222

Ostatni fragment raportu:

Proportion of trace:

LD1 LD2

0.9947 0.0053

pozwala porównywa¢ jako±ci obu funkcji dyskryminacyjnych dzi¦ki tzw. ±ladowi macierzy wariancji i kowariancji (Proportion of trace). Jak wida¢ pierwsza funkcja (LD1) wyja±nia 99.47% wariancji, za± druga jedynie 0.53%. Je±li ana- lizowalismy liniowe funkcje dyskryminacyjne dla dwóch zmiennych obja±niaj¡- cych: Petal.Width oraz Petal.Length to pierwsza funkcja dyskryminacyjna b¦dzie miaªa równanie postaci:

LD1 = -2.402394 * Petal.Width - -1.544371 * Petal.Length a druga:

LD2 =5.042599 * Petal.Width - -2.161222 * Petal.Length

Z kolei wywoªanie komendy > plot(lda.iris,cex=1.0) sprawi, »e uzy- skamy wykres obserwacji z próby ucz¡cej (z podziaªem na klasy) w przestrzeni obu zmiennych kanonicznych (rysunek 6).

7 Zadania do wykonania

1. Prosz¦ utworzy¢ funkcj¦ dyskryminacyjn¡ dla zbioru wine licz¡cego w oryginale 178 obserwacji opisanych 13 atrybutami. Do tworzenia funkcji prosz¦ u»y¢ tylko dwóch zmiennych: V 2 oraz V 8.

2. Równie» dla zbioru wine prosz¦ dokona¢ predykcji przynale»no±ci obser- wacji do klas oraz sporz¡dzi¢ tabele klasykacji.

−6 −4 −2 0 2 4 6

−6−4−20246

LD1

LD2 sssss

s s ss

s s ss ss ss

s ss s

s s

s

s s s

sss s s

s ss

ss s ss s s s s

s s ssss c

c c

c c c

c c

c c

c c

c c

c

cc

c c

c c

c c

c cc c

c c

ccc c c

c c c c

cc c

c c c

c c cc

c

c v

v v

v v

v

v

v v

v v

v v v v v

v

vv v

v v

v v v

v vv v

v v v

v

v

v v

v

v v v v

v

v v

v v

vv v

v

Rysunek 6: Analiza dyskryminacyjna przy dwóch zmiennych kanonicznych

8 Bibliograa

Opracowanie przygotowano w oparciu o prace:

1. J. Koronacki, J. ‚wik: Statystyczne systemy ucz¡ce si¦, wydanie drugie, Exit, Warsaw, 2008. (rozdziaª 3)

2. J. ‚wik, J. Mielniczuk: Statystyczne systemy ucz¡ce si¦ - ¢wiczenia w oparciu o pakiet R, Ocyna Wydawnicza PW, Warszawa, 2009.

3. http://www.ipipan.waw.pl/staff/j.mielniczuk/SSUS-Programy-Dane.

zip

4. http://cran.r-project.org/web/packages/lda/lda.pdf