Optymalny klasyfikator statystyczny

Zadanie klasyfikacji nadzorowanej (por. [Stapor2005], [Stapor2011]) można rozpatrywać w kategoriach probabilistycznych (por. [Greblicki1974], [Berger1985]). Wówczas przestrzeń obserwacji to zbiór wszystkich możliwych wartości wektora cech, natomiast zbiór wszystkich numerów klas \(I=\{1,2,\ldots, L\}\) jest przestrzenią decyzyjną. Prawdopodobieństwa pojawiania się obiektów z poszczególnych klas to prawdopodobieństwa a priori \(P(J=j)=p_j\) dla \(j\in 1,\ldots, L\). Jest to rozkład dyskretnej zmiennej losowej \(J\). Wektor losowy cech \(X\) dla każdego \(j\in\{1,\ldots,k\}\) ma natomiast rozkład prawdopodobieństwa wyrażony gęstością warunkową cech w klasie daną wzorem

\[f(x|j)=f_j(x)\quad\textrm{ dla } x\in\mathbb{R}^d, j\in I\]

Gęstość bezwarunkowa zmiennej \(X\) wyraża się wzorem

\[f(x)=\sum_{j\in J}p_jf_j(x)\]

Przyjmujemy, że \(f(x)>0\) dla każdego \(x\in\mathbb{R}^d\). Obliczmy teraz prawdopodobieństwo warunkowe, że obiekt, któremu odpowiada wektor \(x\) należy do klasy \(j\) tzn.:

\[P(J=j|X=x)\quad\textrm{ dla } j\in\{1,\ldots,L\}, x\in\mathbb{R}\]

Stosując wzór Bayesa otrzymujemy

\[P(J=j|X=x)=\frac{P(X=X|J=j)P(J=j)}{\sum_{j\in\{1,\ldots L\}}P(X=x|J=j)P(J=j)}.\]

Przyjmując oznaczenia

\[\begin{split}p_j=P(J=j)\\ p_j(x)=P(J=j|X=x)\\ f(x)=\sum_{j\in\{1,\ldots L\}}P(X=x|J=j)P(J=j)\end{split}\]

otrzymujemy

\[p_j(x)=\frac{p_jf_j(x)}{f(x)}\]

Prawdopodobieństwo \(p_j(x)\) nazywamy prawdopodobieństwem a posteriori klasy \(j\). Optymalny klasyfikator statystyczny (klasyfikator Bayesa) obiekt reprezentowany przez wektor cech \(x\) przyporządkowuje do tej klasy dla której wartość prawdopodobieństwa a posteriori jest największa. Ponieważ dla dowolnego wektora cech \(x\) wartość \(f(x)\) jest stała a więc reguła decyzyjna klasyfikatora Bayesa jest następująca

\[\phi_{Bayes}(x)=i \quad\textrm{ gdy } p_i(x)f_i(x)=\max_{1\leq j\leq L} p_jf_j(x)\]

Ponieważ logarytm jest funkcją rosnącą a więc powyższą regułę można również zapisać w równoważnej postaci

\[\phi_{Bayes}(x)=i \quad\textrm{ gdy } \ln{p_i(x)f_i(x)}=\max_{1\leq j\leq L} \ln{p_jf_j(x)}\]

Często spotykamy się z przypadkiem kiedy rozkłady wektorów cech pochodzą z wielowymiarowego rozkładu normalnego. Mamy wówczas do czynienia ze szczególnym przypadkiem klasyfikatora Bayesa czyli klasyfikatorem gaussowskim. Funkcje gęstości \(d\)-wymiarowego rozkładu normalnego dane są wzorami

(1)\[f_j(x)=\frac{1}{(2\pi)^{d/2}|\Sigma_j|^{1/2}}exp\left(-\frac{1}{2}(x-m_j)^T\Sigma_j^{-1}(x-m_j)\right)\]

gdzie \(x\) jest \(d\)-wymiarowym wektorem cech, \(m_j\) jest d-wymiarowym wektorem wartości średnich, \(\Sigma_j\) macierzą kowariancji w klasie \(j\). Po podstawieniu (1) uzyskujemy ogólny wzór na funkcje dyskryminacyjne klasyfikatora gaussowskiego

\[g_j(x)=p_j(2\pi)^{-d/2}|\Sigma_j|^{1/2}exp\left(-\frac{1}{2}(x-m_j)^T\Sigma_j^{-1}(x-m_j)\right)\]

Po zlogarytmowaniu otrzymujemy

\[\begin{split}\ln{\left(g_j(x)\right)}&=&\ln{\left(p_j(2\pi)^{-d/2}|\Sigma_j|^{1/2}exp\left(-\frac{1}{2}(x-m_j)^T\Sigma_j^{-1}(x-m_j)\right)\right)}=\\ &=&\ln{p_j}+\ln{\left(2\pi)^{-d/2}\right)}+\ln{\left(|\Sigma_j|^{1/2}\right)}-\frac{1}{2}(x-m_j)^T\Sigma_j^{-1}(x-m_j)=\\ &=&\ln{p_j}-\frac{d}{2}\ln{2\pi}+\frac{1}{2}\ln{|\Sigma_j|}-\frac{1}{2}(x-m_j)^T\Sigma_j^{-1}(x-m_j).\end{split}\]

Pomijając stałą \(\frac{d}{2}\ln{2\pi}\) otrzymujemy funkcje klasyfikujące dla klasyfikatora gaussowskiego

(2)\[g_j(x)=\ln{p_j}+\frac{1}{2}\ln{|\Sigma_j|}-\frac{1}{2}(x-m_j)^T\Sigma_j^{-1}(x-m_j) \quad\textrm{ dla } j\in\{1,\ldots, L\}\]

Z powyższego wzoru wynika, że w ogólnym przypadku funkcja klasyfikująca jest funkcją kwadratową, co z kolei oznacza, że powierzchnie rozdzielające w przestrzeni cech są hiperpłaszczyznami drugiego stopnia a ich typ zależy od parametrów rozkładu cech. W przypadku dwuwymiarowej przestrzeni cech mamy do czynienia z krzywymi płaskimi drugiego stopnia.

Załóżmy, że macierze kowariancji we wszystkich klasach są takie same tzn \(\Sigma_j=\Sigma\) dla każdego \(j\in\{1,\ldots, l\}\). Wówczas po usunięciu jednakowego dla wszystkich klas czynnika \(\frac{1}{2}\ln{|\Sigma|}\) otrzymujemy następującą funkcję dyskryminacyjną

\[g_j(x)=\ln{p_j}-\frac{1}{2}(x-m_j)^T\Sigma_j^{-1}(x-m_j) \quad\textrm{ dla } j\in\{1,\ldots, L\}.\]

Można pokazać, że po dokonaniu w powyższym wzorze szeregu mnożeń a następnie usunięciu stałego dla wszystkich klas czynnika \(x^T\Sigma^{-1}x\) uzyskujemy uproszczoną postać wzoru funkcji dyskryminacyjnych

\[g_j(x)=x^T\Sigma^{-1}m_j-\frac{1}{2}m_j^T\Sigma_j^{-1}m_j+ln(p_j) \quad\textrm{ dla } j\in\{1,\ldots, L\}\]

Wówczas funkcje dyskryminacyjne są funkcjami liniowymi, co z kolei oznacza, że powierzchnie rozdzielające w przestrzeni cech są hiperpłaszczyznami.

Załóżmy dodatkowo, że cechy w poszczególnych klasach są statystycznie niezależne i mają tę samą wariancję \(\sigma^2\) tzn

(3)\[\Sigma_i=\sigma^2I\]

gdzie \(I\) jest macierzą jednostkową stopnia \(L\). Wówczas po podstawieniu (3) do wzoru (2) uzyskujemy

\[g_j(x)=-\frac{1}{2\sigma^2}(x-m_j)^T(x-m_j)+ln(p_j) \quad\textrm{ dla } j\in\{1,\ldots, L\}\]

Stąd po wymnożeniu składnika \((x-m_j)^T(x-m_j)\) i usunięciu wspólnego dla wszystkich klas czynnika \(x^Tx\) otrzymujemy następujący wzór na funkcję dyskryminacyjną

\[g_j(x)=-\frac{1}{2\sigma^2}\left(-2m_j^Tx+m_j^Tm_j\right)+ln(p_j)=w_j^T+w_0\]

gdzie \(w_j=\frac{1}{\sigma^2}m_j\), \(w_0=-\frac{1}{2\sigma^2}m_j^Tm_j+\ln{p_j}\). Również w tym przypadku funkcje dyskryminacyjne są funkcjami liniowymi z czego z kolei wynika, że powierzchnie rozdzielające w przestrzeni cech są hiperpłaszczyznami. Jeśli teraz dodatkowo założymy, że prawdopodobieństwa a priori poszczególnych klas są jednakowe to uzyskamy następująca postać funkcji dyskryminacyjnej

\[g_j(x)=-\frac{1}{2\sigma^2}\left(-2m_j^Tx+m_j^Tm_j\right)\]

Klasyfikator mający funkcje dyskryminacyjne powyższej postaci nazywamy klasyfikatorem minimalno-odległościowym.

Optymalna klasyfikacja bayesowska wymaga znajomości rozkładu zmiennej losowej \(J\) i warunkowej gęstości rozkładu prawdopodobieństwa cech w klasach. W praktycznych zadaniach klasyfikacji tej wiedzy na ogół nie mamy. Spotykamy się z sytuacjami gdzie znana jest postać rozkładu lecz nie są znane jego parametry lub nieznana jest nawet postać rozkładu prawdopodobieństwa. W przypadku gdy znana postać warunkowego rozkładu cech w klasach natomiast nie są znane jego parametry wówczas brakującą wiedzę rekompensujemy informacjami zawartymi w zbiorze uczącym. Zadanie polega na estymacji brakujących parametrów i wstawieniu uzyskanych wartości do funkcji klasyfikujących optymalnego algorytmu Bayesa. Algorytm uzyskany w ten sposób nazywamy Parametrycznym klasyfikatorem Bayesa. Prawdopodobieństwa a priori klas przybliżamy częstościami występowania poszczególnych klas w zbiorze uczącym tzn

\[p_j=\frac{N_j}{N}\]

gdzie \(N\) to liczebność całego zbioru uczącego, natomiast \(N_j\) to liczebność \(j\)-tej klasy. Estymacja parametrów gęstości warunkowych cech w klasach zależy od samej postaci tych gęstości i może odbywać się na przykład metodą największej wiarogodności (por. [Gren1987]).

Przykład 4

Rozważmy zbiór danych Iris. Jest to wielowymiarowy zbiór danych na temat irisów wprowadzony w 1936r. przez Rolanda Fishera. Zbiór ten zawiera \(150\) obiektów z trzech różnych gatunków irysów: setosa, versicolor, virginica. Każdy z obiektów jest charakteryzowany przez cztery cechy: długość (w \(cm\).) listka kielicha kwiatowego, szerokość (w \(cm\)) listka kielicha kwiatowego, długość (w \(cm\)) płatka kwiatu, szerokość (w \(cm\)) płatka kwiatu. Zbiór ten dzielimy, w sposób losowy na dwie części: uczącą i testową w stosunku 70:30. Następnie przeprowadzamy na zbiorze uczącym proces uczenia klasyfikatora Bayesa. Sprawność klasyfikatora to prawie \(98\%\). Jedna próbka ze zbioru testowego została błędnie zaklasyfikowana. Poniżej prezentujemy skrypt w języku python realizujący powyższą procedurę.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from sklearn import datasets
from sklearn.naive_bayes import GaussianNB
from sklearn.cross_validation import train_test_split

#importujemy zbiór danych
iris = datasets.load_iris()

gnb = GaussianNB()

#Zbiór obiektów
X = iris.data

#Wektor poprawnej klasyfikacji obiektów
y = iris.target

y=np.array(y)

#Dzielimy losowo zbiór na dwie części
train, test, train_targets, test_targets = train_test_split(X, y,
                                 test_size=0.30, random_state=42)

#Uczymy klasyfikator
clf = gnb.fit(train, train_targets)

#Testujemy
Z = clf.predict(test)
Następna część - Klasyfikator minimalno-odległościowy