Zagadnienia

7. Klasyfikacja

Zadanie klasyfikacji polega na konstrukcji funkcji (klasyfikatora), która na podstawie zaobserwowanych cech będzie przydzielała obserwację do którejś z wcześniej zdefiniowanych grup. Do estymacji funkcji potrzebne są obserwacje, które już zostały sklasyfikowane, będziemy je nazywać próbą uczącą:

(y_{i},X_{i})\quad i\in\{ 1,\ldots,n\}\text{ niezależne obserwacje, }X_{i}\in\mathbb{R}^{p}\quad y_{i}\in\{ 1,\ldots,K\}.

Dane X_{i} oznaczają zaobserwowane cechy, y_{i} grupę, do której obserwacja została zaklasyfikowana.

C^{j} oznaczaja zbiór tych indeksów i, że y_{i}=j; C^{1},\ldots,C^{K} są rozłączne o licznościach odpowiednio n_{1},\ldots,n_{K}.

Funkcja wiarygodności dla opisanych danych wyraża się wzorem:

L(\pi _{1},\ldots,\pi _{K},\theta _{1},\ldots,\theta _{K})=\prod _{{i=1}}^{n}g(y_{i},X_{i})=\prod _{{i=1}}^{n}\pi _{{y_{i}}}f_{{\theta _{{y_{i}}}}}(X_{i})
\text{pod warunkiem: }\sum _{{k=1}}^{K}\pi _{k}=1.

Logwiarygodność to logarytm funkcji wiarygodności:

l=\log L;
l=\sum _{{i=1}}^{n}\log\pi _{{y_{i}}}+\sum _{{i=1}}^{n}\log f_{{\theta _{{y_{i}}}}}(X_{i}).

Zadanie klasyfikacji dzielimy na dwa kroki:

  1. Estymujemy parametry \pi _{1},\ldots,\pi _{K} oraz \theta _{1},\ldots,\theta _{K} na podstawie zaobserwowanych par (y_{i},X_{i}) przy użyciu metody największej wiarygodności. Parametry \pi _{k} możemy interpretować jako prawdopodobieństwa przynależności do danej grupy danych, a \theta _{k} jako parametry rozkładu w danej grupie (na przykład dla wielowymiarowego rozkładu normalnego, byłyby to średnia \mu i macierz kowariancji \Sigma).

  2. Obserwujemy nowe cechy X_{{n+1}} i przyporządkowujemy im \widehat{y}_{{i+1}} na podstawie zbudowanego przez nas klasyfikatora. Będziemy go także nazywać regułą decyzyjną.

Maksymalizujemy funkcję wiarygodności pod warunkiem \sum _{{k=1}}^{K}\pi _{k}=1 przy użyciu metody mnożników Lagrange'a:

\max _{{\pi _{1},\ldots,\pi _{K},\theta _{1},\ldots,\theta _{K}}}F(\pi _{1},\ldots,\pi _{K},\theta _{1},\ldots,\theta _{K})=
\max _{{\pi _{1},\ldots,\pi _{K},\theta _{1},\ldots,\theta _{K}}}n_{1}\log\pi _{1}+\ldots+n_{K}\log\pi _{K}-\lambda(\sum _{{k=1}}^{K}\pi _{k}-1)+
+\sum _{{i\in C^{1}}}\log f_{{\theta _{1}}}(X_{i})+\ldots\sum _{{i\in C^{K}}}\log f_{{\theta _{K}}}(X_{i}).

Liczymy estymatory \widehat{\pi}_{1},\ldots,\widehat{\pi}_{K}:

\frac{\partial F}{\partial\pi _{k}}=\frac{n_{k}}{\pi _{k}}-\lambda=0\quad\forall k=1\ldots,K; (7.1)
\frac{\partial F}{\partial\lambda}=\sum _{{k=1}}^{K}\pi _{k}-1=0; (7.2)

Z równań 7.1 otrzymujemy:

\forall k=1\ldots,K\quad\frac{n_{k}}{\lambda}=\pi _{k}.

Sumujemy po k korzystając z równania 7.2 :

\frac{n}{\lambda}=1\quad\Rightarrow\quad\lambda=n\quad\Rightarrow\quad\widehat{\pi}_{k}=\frac{n_{k}}{n}.

Estymację parametrów \theta _{k} odłożymy do dalszej części wykładu.

7.1. Optymalna reguła decyzyjna

Zobaczmy teraz, jak można zdefiniować optymalny klasyfikator w zależności od funkcji straty karzącej za błędne sklasyfikowanie danych.

Definicja 7.1

Funkcja straty to funkcja przyporządkowująca nieujemną wielkość kary poprzez porównanie prawdy (założymy chwilowo, że ją znamy) do podjętej decyzji (wyliczonego estymatora):

L:\underbrace{K}_{{prawda}}\times\underbrace{K}_{{decyzja}}\rightarrow\underbrace{\mathbb{R}^{+}}_{{kara}}.

Przykładową funkcją kary dla ciągłego y jest L(y,\widehat{y})=(y-\widehat{y})^{2}.

Mając wyestymowaną regułę decyzyjną sensownym jest rozpatrywanie średniej straty dla naszego klasyfikatora:

Definicja 7.2

Ryzyko reguły decyzyjnej dla d:\mathbb{R}^{p}\rightarrow K:

\text{Ryzyko}=\text{średnia strata reguły decyzyjnej }d=
=\text{R}(d)=\sum _{{y\in K}}\int _{{\mathbb{R}^{p}}}L(y,d(X))g(y,X)dX=

gdzie g(y,X) jest gęstością łącznego rozkładu danych. Z twierdzenia o prawdopodobieństwie całkowitym:

=\sum _{{k=1}}^{K}\left[\int _{{\mathbb{R}^{p}}}L(k,d(X))f(X|k)dX\right]\pi _{k}.
Definicja 7.3

Optymalna reguła decyzyjna d_{*} to taka reguła decyzyjna, że

\forall d\quad\text{R}(d_{*})\leq\text{R}(d).
Definicja 7.4

Reguła bayesowska d_{B}(X) to reguła decyzyjna, która lokalnie dla danego X spełnia warunek:

d_{B}(X)=\text{argmin}_{{1\leq l\leq K}}\mathbb{E}_{{y|X}}L(y,l)=
=\text{argmin}_{l}\sum _{{k=1}}^{K}L(k,l)p(k|X)=

ze wzoru Bayesa:

=\text{argmin}_{l}\left[\sum _{{k=1}}^{K}L(k,l)\frac{\pi _{k}f(X|k)}{\sum _{{s=1}}^{K}\pi _{s}f(X|s)}\right]=
=\text{argmin}_{l}\left[\sum _{{k=1}}^{K}L(k,l)\pi _{k}f(X|k)\right].
Stwierdzenie 7.1
\text{R}(d_{B})=\text{R}(d_{*}).

Reguła bayesowska jest optymalną regułą decyzyjną.

Dla dowolnej reguły decyzyjnej d zachodzi:

\text{R}(d)=\sum _{{k=1}}^{K}\left[\int _{{\mathbb{R}^{p}}}L(k,d(X))f(X|k)dX\right]\pi _{k}=
=\int _{{\mathbb{R}^{p}}}\left[\sum _{{k=1}}^{K}L(k,d(X))\pi _{k}f(X|k)\right]dX\geq
\geq\int _{{\mathbb{R}^{p}}}\left[\min _{{1\leq l\leq k}}\sum _{{k=1}}^{K}L(k,l)\pi _{k}f(X|k)\right]dX=
=\int _{{\mathbb{R}^{p}}}\left[\sum _{{k=1}}^{K}L(k,d_{B}(X))\pi _{k}f(X|k)\right]dX=\text{R}(d_{B}).

7.2. Wielowymiarowy rozkład normalny

W dalszej części wykładu będziemy zakładać, że f_{{\theta _{i}}} mają rozkłady normalne. Dlatego przyjrzyjmy się bliżej własnościom wielowymiarowego rozkładu normalnego i estymacji jego parametrów metodą najwiękzej wiarygodności.

Definicja 7.5

Wektor losowy X=(x_{1},\ldots,x_{p}) ma rozkład wielowymiarowy normalny w \mathbb{R}^{p} jeśli \forall u\in\mathbb{R}^{p} u^{T}X ma rozkład normalny w \mathbb{R}. Oznaczmy ten rozkład poprzez \mathcal{N}(\mu,\Sigma), gdzie \mu=\mathbb{E}x, \Sigma=\text{Var}(X).

Twierdzenie 7.1

Jeżeli X ma rokład normalny w \mathbb{R}^{p}, to \forall a\in\mathbb{R}^{k} i macierzy A wymiaru k\times p, AX+a ma rozkład normalny w \mathbb{R}^{k}.

\forall\text{ }u\in\mathbb{R}^{k}\quad u^{T}(AX+a)=(u^{T}a)X+u^{T}a.
Wniosek 7.1

Rozkłady brzegowe wielowymiarowego rozkładu normalnego są normalne w odpowiednich podprzestrzeniach \mathbb{R}^{p}.

Twierdzenie 7.2

Fcja charakterystyczna zmiennej losowej X o rozkładzie normalnym w \mathbb{R}^{p} jest postaci:

\varphi _{X}(t)=e^{{it^{T}\mu-\frac{1}{2}t^{T}\Sigma t}}. (7.3)

Także na odwrót: jeżeli \Sigma jest symetryczną macierzą dodatnio określoną o wymiarach p\times p, to \varphi _{X} określona w równaniu 7.3 jest funkcją charakterystyczną wektora losowego o rozkładzie normalnym w \mathbb{R}^{p}.

Wniosek 7.2

Dowolna macierz symetryczna dodatnio określona o wymiarach p\times p jest macierzą kowariancji wektora losowego o rokładzie normalnym w \mathbb{R}^{p}.

Twierdzenie 7.3

Gęstość wielowymiarowego rozkładu normalnego \mathcal{N}(\mu,\Sigma):

f(X)=\frac{1}{(2\pi)^{{\frac{p}{2}}}(\underbrace{\text{det}(\Sigma)}_{{=|\Sigma|}})^{{\frac{1}{2}}}}\exp\left[-\frac{1}{2}(X-\mu)^{T}\Sigma^{{-1}}(X-\mu)\right].
Twierdzenie 7.4

Jeżeli X ma rozkład normalny w \mathbb{R}^{p}: \mathcal{N}(\mu,\Sigma), to współrzędne wektora X są niezależne \Leftrightarrow \Sigma jest diagonalna. Dla rozkładu normalnego brak korelacji oznacza niezależność.

Twierdzenie 7.5

Jeżeli X\sim\mathcal{N}(\mu,\sigma^{2}\mathbb{I}), C jest macierzą ortonormalną o wymiarach p\times p, to:

CX\sim\mathcal{N}(C\mu,C\sigma^{2}\mathbb{I}C^{T})=\mathcal{N}(C\mu,\sigma^{2}\underbrace{CC^{T}}_{{=\mathbb{I}}})=\mathcal{N}(C\mu,\sigma^{2}\mathbb{I}).

7.2.1. Estymatory największej wiarygodności dla rozkładu normalnego \mathcal{N}(\mu,\Sigma)

Niech X_{1},\ldots,X_{n} będą niezależnymi wektorami losowymi z p-wymiarowego rozkładu \mathcal{N}(\mu,\Sigma). Znajdziemy estymatory dla parametrów \mu i \Sigma. Łączna funkcja wiarygodności dla n wektorów losowych:

L(\mu,\Sigma)=\prod _{{i=1}}^{n}\frac{1}{(2\pi)^{{\frac{p}{2}}}|\Sigma|^{{\frac{1}{2}}}}\exp\left[-\frac{1}{2}(X_{i}-\mu)^{T}\Sigma^{{-1}}(X_{i}-\mu)\right].

Najpierw szukamy estymatora \widehat{\mu}; w tym celu opuszczamy wszystkie wyrazy nie zalezeżące od \mu, które by się wyzerowały po policzeniu pochodnej. Dla prostoty obliczeń maksymalizujemy podwojoną logwiarygodność:

2\log(L(\mu))=2l(\mu)=n\mu^{T}\Sigma^{{-1}}-2n\mu^{T}\Sigma^{{-1}}\overline{X},

gdzie \overline{X}=\frac{1}{n}\sum _{{i=1}}^{n}X_{i}.

Przypomnijmy fakt:

Lemat 7.1

Oznaczmy: a,b – wektory tej samej długości p, A macierz o wymiarach p\times p.

\frac{\partial(a^{T}b)}{\partial a}=\frac{\partial(b^{T}a)}{\partial a}=b.
\frac{\partial(b^{T}Ab)}{\partial b}=\underbrace{(A+A^{T})b}_{{=2Ab\text{ jeśli $A$ symetryczna}}}.

Skorzystajmy z lematu 7.1 żeby obliczyć pochodną logwiarygodności:

\frac{1}{n}\frac{\partial(2l(\mu))}{\partial\mu}=2\Sigma^{{-1}}\mu-2\Sigma^{{-1}}\overline{X}=0,

stąd

\widehat{\mu}=\overline{X},

czyli estymatorem największej wiarygodności dla średniej rozkładu normalnego jest średnia arytmetyczna obserwacji.

Ponieważ optymalne \widehat{\mu} nie zależy od \Sigma, przy obliczaniu \widehat{\Sigma} możemy wstawić \overline{X} za \mu. Maksymalizujemy po \Sigma wyrażenie:

L(\overline{X},\Sigma)\propto|\Sigma|^{{-\frac{n}{2}}}\exp\left[-\frac{1}{2}\sum _{{i=1}}^{n}(X_{i}-\overline{X})^{T}\Sigma^{{-1}}(X_{i}-\overline{X}).\right]

Symbol \propto oznacza proporcjonalność, możemy opuścić wszystkie stałe, które nie wpływają na wynik optymalizacji.

Ponieważ (X_{i}-\overline{X})^{T}\Sigma^{{-1}}(X_{i}-\overline{X}) jest liczbą, a \text{tr}(\text{liczba})=\text{liczba}, oraz \text{tr}(AB)=\text{tr}(BA), otrzymujemy:

L(\overline{X},\Sigma)\propto|\Sigma|^{{-\frac{n}{2}}}\exp\left[-\frac{1}{2}\sum _{{i=1}}^{n}\text{tr}\{(X_{i}-\overline{X})(X_{i}-\overline{X})^{T}\Sigma^{{-1}}\}\right]=

ślad macierzy jest funkcją liniową argumentu, więc zachodzi:

=|\Sigma|^{{-\frac{n}{2}}}\exp\left[-\frac{1}{2}\text{tr}\{\underbrace{(X_{i}-\overline{X})(X_{i}-\overline{X})^{T}}_{{=S}}\Sigma^{{-1}}\}\right]=

pomnóżmy i podzielmy przez |S|^{{\frac{n}{2}}}

=|S|^{{-\frac{n}{2}}}|S\Sigma^{{-1}}|^{{\frac{n}{2}}}\exp\left[-\frac{1}{2}\text{tr}\{ S\Sigma^{{-1}}\}\right].

Ponieważ |S|^{{-\frac{n}{2}}} nie zależy od \Sigma, możemy to wyrażenie opuścić. Podstawmy B=S\Sigma^{{-1}}:

L(\overline{X},B)\propto|B|^{{\frac{n}{2}}}\exp\left[-\frac{1}{2}\text{tr}B\right].
Lemat 7.2

Dla macierzy kwadratowej A o wymiarach p\times p zachodzi:

\text{det}(A)=\prod _{{i=1}}^{p}\lambda _{i},
\text{tr}(A)=\sum _{{i=1}}^{p}\lambda _{i},

gdzie \lambda _{i} to wartości własne macierzy.

Korzystając z lematu 7.2:

|B|^{{\frac{n}{2}}}\exp\left[-\frac{1}{2}\text{tr}B\right]=\prod _{{j=1}}^{p}\lambda _{j}^{{\frac{n}{2}}}e^{{-\frac{1}{2}\lambda _{j}}}

Zmaksymalizujmy to wyrażenie po każdej wartości własnej \lambda _{j}, co sprowadza się do maksymalizacji po \lambda funkcji:

F(\lambda)=\lambda^{{\frac{n}{2}}}e^{{-\frac{1}{2}\lambda}};
\log F(\lambda)=\frac{n}{2}\log\lambda-\frac{1}{2}\lambda;
\frac{\partial(\log F(\lambda))}{\partial\lambda}=\frac{n}{2\lambda}-\frac{1}{2}=0;

skąd \widehat{\lambda}=\widehat{\lambda _{1}}=\ldots=\widehat{\lambda _{p}}=n.

Macierzą o wszystkich wartościach własnych równych n jest n\mathbb{I}:

B=S\Sigma^{{-1}}=n\mathbb{I},

skąd:

\widehat{\Sigma}=\frac{1}{n}S=\frac{1}{n}(X_{i}-\overline{X})(X_{i}-\overline{X})^{T},

czyli estymatorem największej wiarygodności dla macierzy kowariancji rozkładu normalnego jest obciążony estymator próbkowy macierzy kowariancji.

7.3. Klasyfikacja w modelu normalnym

Zrobimy dwa założenia dotyczące rozważanego wcześniej klasyfikatora:

  1. Funkcja straty jest postaci:

    L(k,l)=\textbf{1}_{{k\neq l}}.
  2. W każdej z grup dane pochodzą z rozkładu normalnego, czyli f_{{\theta _{k}}} to gęstość rozkładu normalnego, \theta _{k}=(\mu _{k},\Sigma _{k}).

Dal zadanej funkcji straty optymalna (bayesowska) reguła decyzyjna będzie miała postać:

d_{B}(X)=\text{argmin}_{l}\left[\sum _{{k=1}}^{K}L(k,l)\pi _{k}f(X|k)\right]=\text{argmin}_{l}\left[\sum _{{k=1}}^{K}\textbf{1}_{{k\neq l}}\pi _{k}f(X|k)\right]=
=\text{argmin}_{l}\left[\underbrace{\sum _{{k=1}}^{K}\pi _{k}f(X|k)}_{{\text{nie zależy od wyboru }l}}-\pi _{l}f(X|l)\right]=\text{argmax}_{l}\left[\pi _{l}f(X|l)\right].

Znamy już postać szukanego klasyfikatora, potrzebujemy jeszcze estymatorów dla występujących w nim parametrów. Wiemy jak wyglądają estymatory \widehat{\pi}_{k}:

\widehat{\pi}_{k}=\frac{n_{k}}{n},\quad n_{k}=\sum _{{i=1}}^{n}\textbf{1}(y_{i}=k).

Estymatory największej wiarygodności dla parametrów \theta _{k} przy założeniu normalności rozkładów w grupach są postaci:

\widehat{\mu _{k}}=\frac{1}{n_{k}}\sum _{{i\in C^{k}}}X_{i}=\frac{\sum _{{i=1}}^{n}X_{i}\textbf{1}(y_{i}=k)}{\sum _{{i=1}}^{n}\textbf{1}(y_{i}=k)};
\widehat{\Sigma _{k}}=\frac{1}{n_{k}}\sum _{{i\in C^{k}}}(X_{i}-\overline{X^{k}})(X_{i}-\overline{X^{k}})^{T},

gdzie \overline{X^{k}} oznacza wektor średnich obserwacji dla X_{i}\in C^{k}.

Dla X niezależnego od próby uczącej: (y_{1},X_{1}),\ldots,(y_{n},X_{n}) estymator reguły decyzyjnej ma postać:

\widehat{d}(X)=\text{argmax}_{{1\leq l\leq K}}\widehat{\pi _{l}}\left[f_{{\widehat{\mu}_{l},\widehat{\Sigma}_{l}}}(X)\right].

7.3.1. Kwadratowa (qda) i liniowa (lda) funkcja klasyfikacyjna

W zależności od założeń dotyczących parametrów, możemy otrzymać klasyfikator będący różną funkcją swojego argumentu X: albo kwadratową albo liniową.

Kwadratowa funkcja klasyfikacyjna (qda) nie wymaga dodatkowych założeń o parametrach:

d(X)=\text{argmax}_{l}\left[\pi _{l}f_{{\mu _{l},\Sigma _{l}}}\right]=
=\text{argmax}_{l}\left[\frac{\pi _{l}}{(2\pi)^{{\frac{p}{2}}}|\Sigma _{l}|^{{\frac{1}{2}}}}\exp\left\{-\frac{1}{2}(X-\mu _{l})^{T}\Sigma _{l}^{{-1}}(X-\mu _{l})\right\}\right]=

po opuszczeniu wyrażeń niezależnych od l i zlogarytmowaniu:

=\text{argmax}_{l}\left[\log(\pi _{l})-\frac{1}{2}\log(|\Sigma _{l}|)-\frac{1}{2}(X-\mu _{l})^{T}\Sigma _{l}^{{-1}}(X-\mu _{l})\right],

czyli kwadratowa funkcja argumentu X.

Liniowa funkcja klasyfikacyjna (lda) wymaga założenia:

\Sigma _{1}=\ldots=\Sigma _{K}=\Sigma.

Dzięki niemu mamy podwójny zysk obliczeniowy: o K-1 parametrów mniej do wyestymowania i liniową funkcję optymalizowaną:

d(X)=\text{argmax}_{l}\left[\log(\pi _{l})-\frac{1}{2}\log(|\Sigma _{l}|)-\frac{1}{2}(X-\mu _{l})^{T}\Sigma _{l}^{{-1}}(X-\mu _{l})\right]=

ponieważ \log(|\Sigma|) oraz X^{T}\Sigma^{{-1}}X nie zależy od l,

=\text{argmax}_{l}\left[\log(\pi _{l})+X^{T}\Sigma^{{-1}}\mu _{l}-\frac{1}{2}\mu _{l}^{T}\Sigma^{{-1}}\mu _{l}\right].

7.4. Metody porównywania klasyfikatorów

Chcemy znaleźć taką metodę porónywania, żeby każdą obserwację spośród (y_{1},X_{1}),\ldots,(y_{n},X_{n}) wykorzystać do uczenia i testu, ale tak żeby testować tylko na tych obserwacjach, które nie były brane pod uwagę przy uczeniu klasyfikatorów.

Kroswalidacja m-krotna (walidacja krzyżowa) polega na podziale danych na m części (popularnymi wyborami są m=5, m=10): m-1 będzie tworzyć próbę uczącą, ostatnia będzie próbą testową. Estymujemy klasyfikatory na próbie uczącej, porównujemy metody na próbie testowej. Powtarzamy procedurę m razy tak, żeby każda z części była próbą testową. Dokładniej:

  1. Permutujemy obserwacje. Jeżeli dane mają jakąś strukturę, na przykład można je podzielić na klasy, permutujemy obserwacje w klasach.

  2. Dzielimy próbę na m części tak, żeby w każdej z grup było po tyle samo obserwacji z każdej klasy.

  3. Uczymy klasyfikatory na próbie uczącej - estymujemy parametry.

  4. Porównujemy metody na próbie testowej (np. poprzez estymację prawdopodobieństwa poprawnej predykcji)

Definicja 7.6

Prawdopodobieństwo poprawnej predykcji to dla danego klasyfikatora \mathbb{P}(d(X)=y). Np. jeżeli funkcja straty wyraża się wzorem L(k,l)=\textbf{1}_{{k\neq l}}, możemy estymować prawdopodobieństwo poprawnej predykcji dla konkretnej próby treningowej i testowej następująco:

\widehat{ppp}_{i}=\frac{\sum _{{i\in\text{próba testowa}}}{\textbf{1}(d(X_{i})=y_{i})}}{\sum _{{i\in\text{próba testowa}}}1},

gdzie d jest klasyfikatorem wyestymowanym na podstawie próby uczącej. Uśrednione \widehat{ppp} jest dobrą metodą porównywania klasyfikatorów:

\widehat{ppp}=\frac{\sum _{{i=1}}^{m}\widehat{ppp}_{i}}{m}.

7.5. Przykłady w programie R

Klasyfikacja:

Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.

Projekt współfinansowany przez Unię Europejską w ramach Europejskiego Funduszu Społecznego.

Projekt współfinansowany przez Ministerstwo Nauki i Szkolnictwa Wyższego i przez Uniwersytet Warszawski.