Statystyczna Eksploracja Danych

LABORATORIUM 2

GĘSTOŚĆ ROZKŁADÓW 2D

Głownym zdaniem ćwiczeń będzie dokonanie dyskryminacji za pomocą metod LDA oraz QDA, które zostały w przedstawione na Wykładzie 2. Zostaną do tego wykorzystane dedykowane funkcje z bibliotek MASS oraz klaR, jednak aby w pełni zaznajomić się z tymi metodami warto będzie wyreślić najprostsze hiperpłaszczyzny dyskryminujące "na piechotę". W tym celu należy poznać sposoby wyświetlania poziomic rozkładów 2D za pomocą funkcji contour(), a co za tym idzie również rozpinania siatki na zakresach współrzędnych.

Zaczynamy od uruchomienia dwóch bibliotek: MASS oraz mvtnorm (tę drugą należy zainstalować). Pierwszą już znamy (m.in. funkcje ginv() oraz mvrnorm()), natomiast drugą wykorzystamy do wyznaczania gęstości prawdopodobieństwa wielomwymarowego rozkładu Gaussa (funkcja dmvnorm()).

library(MASS)
library(mvtnorm)

Podobnie jak podczas poprzedniego laboratorium generujemy punkty z dwóch rozkładów Gaussa o zadanych wektorach \(\mathbf{m}_1\) i \(\mathbf{m}_2\) oraz macierzach \(\mathbf{S}_1\) i \(\mathbf{S}_2\)

# Macierze kowariancji dla klas
S1 <- matrix(c(4,0,0,4),2,2)
S2 <- matrix(c(2,0,0,2),2,2)

# Wartości oczekiwane
mt1 <- c(-3, -1)
mt2 <- c(2, 2)

# Liczba punktów w klasach
n1 <- 40
n2 <- 30
n <- n1 + n2

# Generowanie rozkładów
X1 <- mvrnorm(n1, mt1, S1)
X2 <- mvrnorm(n2, mt2, S2)

# Zamiana na ramki danych
X1 <- data.frame(X1); colnames(X1) <- c("x", "y")
X2 <- data.frame(X2); colnames(X2) <- c("x", "y")

Następnie wykreślamy punkty

# Nanoszenie punktów na wykres
plot(X1$x, X1$y, ylim = c(-8,8), xlim = c(-8,8), xlab = "X", ylab = "Y", pch = 19, col = "blue", font = 2)
abline(v = 0, h = 0, col = "gray")
points(X2$x, X2$y, pch = 19, col = "orange")

otrzymując efekt podobny do poniższego

Rysunek 2.1

Tworzenie dwuwymiarowego rozkładu gęstości prawdopodobieństwa (PDF) jest oczywiście ideowo tożsame z przypadkiem jednowymiarowym: korzystając z odpowiednich funkcji (dnorm() dla 1D oraz dmvnorm() dla 2D), do których przekazujemy parametry rozkładu (\(\mu, \sigma\) dla 1D, \(\mathbf{m},\mathbf{S}\) dla 2D) otrzymujemy gęstość prawdopodobieństwa w danym punkcie. Różnica polega na tym, że w przypadku 2D musimy przekazać do funkcji obie współprzędne \((x_i,y_i)\)każdego punktu \(\mathbf{x}_i\). Rzecz jasna, można po prostu zadeklarować macierz o odpowiednich rozmiarach wypełnioną zerami, a następnie korzystając z dwóch pętli wypełniać ją, jednak takie podejście jest mało "eleganckie". Zamiast tego:

  1. zdefiniujemy zakresy współprzędnych,
  2. rozepniemy na nich siatkę (funkcja expand.grid()),
  3. wyznaczymy wartości PDF 2D Gaussa dla każdego z tych punktów i zapiszemy je do macierzy.

Poniższy kod realizuje te procedury.

# Definiowanie zakresów współprzędnych
xp <- seq(-10, 10, 0.1)
yp <- seq(-10, 10, 0.1)

# Rozpięcie siatki na współprzędnych
gr <- expand.grid(x = xp, y = yp)

# Wyznaczenie gęstości prawdopodobieństwa dla każdego punktu siatki
# i zapisanie efektów do macierzy
X1.pdf <- matrix(dmvnorm(gr, mt1, S1), length(xp))
X2.pdf <- matrix(dmvnorm(gr, mt2, S2), length(xp))

Ostatnim punktem jest wykreślenie warstwic rozkładu za pomocą funkcji contour() dla obu klas

# Wykreślenie warstwic
contour(xp, yp, X1.pdf, add = TRUE, col = "blue")
contour(xp, yp, X2.pdf, add = TRUE, col = "orange")

otrzymując efekt podobny do poniższego

Rysunek 2.2

METODA LDA

Przechodzimy teraz do liniowej analizy dyskryminacji (LDA). Będziemy korzystać z funkcji lda() lecz wcześniej musimy umieścić wszystkie dane w jednej ramce danych, jak również przypisać im przynależność do klas.

# Przypisanie klas
X1$class <- 1; X2$class <- 2

# "Sklejenie" danych do jednej ramki
# oraz przekształcenie typu zmiennej przechowującej klasy
data <- rbind(X1, X2); data$class <- factor(data$class)

Teraz możemy wykonać już funkcję lda() i zapisać efekt jej działania do zmiennej. Sama funkcja opowiada jedynie za proces konstrukcji klasyfikatora, natomiast do przewidywania przynależności dla nowych danych korzystamy z funkcji predict() (funkcja predict() jest podobnie przeciążana jak np. plot() - będziemy korzystać z niej również dla innych klasyfikatorów).

# Metoda LDA
data.lda <- lda(class ~ x + y, data)

# Przewidywanie klas za pomocą metody LDA
# korzystamy z tych samych danych co przy uczeniu
# czyli powtórne podstawienie
data.lda.pred <- predict(data.lda, data)

Łatwo można porównać wyniki klasyfikacji za pomocą funkcji table(). W ten sposób tworzymy właśnie macierz pomyłek, musimy jedynie przekazać zarówno oryginalne klasy (ze zmiennej data) jak również jedno z pól zmiennej data.lda.pred (pole class).

# Budowanie macierzy pomyłek
CM <- table(data$class, data.lda.pred$class)

Kluczowa informacja, dotycząca przynależności do klas jest zawarta w polu $posterior. Znajdują się tam prawdopodobieństwa a posteriori przynależności do klas dla wszystkich obserwacji. Klasa, dla której prawdopodobieństwo jest największe jest ostatecznie przypisywana do obserwacji i umieszczana w polu $class.

Spróbujmy teraz odtworzyć wartości tych prawdopodobieństw; w przypadku dwóch klas de facto wystarczy nam informacja o \(p(1|\mathbf{x})\), które jest dane wzorem (reguła Bayesa)

\(p(1|\mathbf{x}) = \frac{\pi_1 p(\mathbf{x} | 1)}{p(\mathbf{x})} = \frac{\pi_1 p(\mathbf{x} | 1)}{\pi_1 p(\mathbf{x} | 1) + \pi_2 p(\mathbf{x} | 2)}\)

Prawdopodobieństwa a priori \(\pi_1\) oraz \(\pi_2\) estymujemy jako liczności klas podzielone przez całkowitą liczbę obserwacji, natomiast gęstości prawdopodobieństwa \(p(\mathbf{x} | 1)\) i \(p(\mathbf{x} | 2)\) są dwuwymiarowymi rozkładami Gaussa o średnich \(\mathbf{m}_1\) i \(\mathbf{m}_2\) estymowanych z danych. Istotną założeniem metody LDA jest to, że macierze kowariancji w obu klasach powinny być takie same. W praktyce uzyskuje się to poprzez uśrednienie wyestymowanych macierzy kowariancji \(\mathbf{S}_1\) i \(\mathbf{S}_2\)

\(\mathbf{S} = \frac{1}{n - 2} \sum\limits_{k=1}^{2}(n_k - 1)\mathbf{S}_k\)

Po tych wszystkich uwagach możemy teraz wyznaczyć wartości prawdopodobieństwa a posteriori przynależności do 1 klasy i porównać je z tymi, które otrzymaliśmy bezpośrednio z funkcji predict()

# Funkcja do wyznaczania prawdopodobieństw a posteriori w metodzie LDA
f.lda <- function(X, m1, m2, S, pi1, pi2) {
return(pi1 * dmvnorm(X, m1, S) / (pi1 * dmvnorm(X, m1, S) + pi2 * dmvnorm(X, m2, S)))
}

# Estymowane wartości oczekiwane
me1 <- apply(X1[,1:2], 2, mean)
me2 <- apply(X2[,1:2], 2, mean)

# Estymowane macierze kowariancji
Se1 <- cov(X1[,1:2])
Se2 <- cov(X2[,1:2])

# Macierz kowariancji do metody LDA
# tożsama z macierzą W
Se <- ((n1 - 1) * Se1 + (n2 - 1) * Se2) / (n1 + n2 - 2)

# Prawdopodbieństwa a priori
pi1 <- n1 / (n1 + n2)
pi2 <- n2 / (n1 + n2)

# Porównanie prawdopodobieństw a posteriori
data.comp <- cbind(f.lda(data[,1:2], me1, me2, Se, pi1, pi2), data.lda.pred$posterior[,1])

Suche liczby nigdy nie wytrzymują porównania z wykresem. Biblioteka klaR (do zainstalowania) oferuje funkcję partimat(), za pomocą której można wykreślić rozgraniczenie klas (prostą dyskryminującą). Niestety, ma ona swoje mankamenty i zamiast niej użyjemy funkcji drawparti() (w zasadzie obie powinny dawać to samo), ułatwiając sobie życie za pomocą komendy with().

library(klaR)

# Rysowanie prostej rozdzielającej, punktów etc
# Dla większej dokładności użyć opcji prec = 200
with(data, drawparti(class, x, y, xlab = "X", ylab = "Y", font = 2))

Rysunek 2.3

Na koniec wykorzystamy obliczone przez nas prawdopodobieństwa a posteriori: prosta rozgraniczająca klasy powinna być po prostu warstwicą 0.5 prawdopodobieństwa a posteriori \(p(1|\mathbf{x})\).

# Porównanie z wartościami otrzymanymi z rozkładów
contour(xp, yp, matrix(f.lda(gr, me1, me2, Se, pi1, pi2), length(xp)), add = T, levels = 0.5, lwd = 2, lty = 2, col = "blue")

Rysunek 2.4

Warto pamiętać, że to samo powinniśmy otrzymać prosto z teorii, przekształcając równanie prostej dyskryminującej

\(\ln \frac{\pi_1}{\pi_2} - \frac{1}{2} (\mathbf{m}_1-\mathbf{m}_2)^{T} \mathbf{S}^{-1} (\mathbf{m}_1+ \mathbf{m}_2) + (\mathbf{m}_1-\mathbf{m}_2)^{T} \mathbf{S}^{-1} \mathbf{x} = 0\)

do wygodniejszej postaci

\(a_x x + a_y y + b = 0\), gdzie
\(b = \ln \frac{\pi_1}{\pi_2} - \frac{1}{2}\mathbf{a}^{T}(\mathbf{m}_1 + \mathbf{m}_2)\)
\(\mathbf{a}^{T} = (\mathbf{m}_1-\mathbf{m}_2)^{T} \mathbf{S}^{-1}\)

czyli

# Wyznaczenie wektora a oraz wyrazu b
a <- t(me1 - me2) %*% ginv(Se)
b <- log(pi1 / pi2) - 0.5 * a %*% (me1 + me2)

# Wyznaczenie wyrazu wolnego i współczynnika kierunkowego
B <- -b / a[2]
A <- -a[1] / a[2]

# Wykreślanie prostej
abline(B, A, col = "red", lwd = 2)

co skutkuje

Rysunek 2.5

potwierdzając, że wszystkie trzy podejścia są tożsame.

METODA QDA

W przypadku kwadratowej analizy dyskryminacji (QDA) korzystamy z funkcji qda(), natomiast dalsze instrukcje są analogiczne do poprzednich

# Metoda QDA
data.qda <- qda(class ~ x + y, data)

# Przewidywanie klas za pomocą metody QDA
data.qda.pred <- predict(data.qda, data)

# Budowanie macierzy pomyłek
CM <- table(data$class, data.qda.pred$class)

Do wykreślenia granic klas funkcja drawparti() musi zostać wyposażona w opcję method="qda"

# Rysowanie prostej rozdzielającej, punktów etc
with(data, drawparti(class, x, y, method = "qda", xlab = "X", ylab = "Y", font = 2))

Rysunek 2.6

Aby znów wykorzystać podejście bezpośredniego wyznaczania prawdopodobieństwa a posteriori, należy pamiętać, że w metodzie QDA zakłada się różne macierze kowariancji \(\mathbf{S}_1\) i \(\mathbf{S}_2\), stąd musimy nieco zmodyfikować oryginalną funkcję

# Funkcja do wyznaczania prawdopodobieństw a posteriori w metodzie QDA
f.qda <- function(X, m1, m2, S1, S2, pi1, pi2) {
return(pi1 * dmvnorm(X, m1, S1) / (pi1 * dmvnorm(X, m1, S1) + pi2 * dmvnorm(X, m2, S2)))
}

# Porównanie z wartościami otrzymanymi z rozkładów
contour(xp, yp, matrix(f.qda(gr, me1, me2, Se1, Se2, pi1, pi2), length(xp)), add = T, levels = 0.5, lwd = 2, lty = 2, col = "blue")

Rysunek 2.7

NAIWNY BAYES

Warto wspomnieć jaka jest (bardzo istotna!) różnica pomiędzy klasyfikatorem Bayesa, który stanowi podstawę metod LDA i QDA oraz naiwnym klasyfikatorem Bayesowskim. Otóż w tym ostatnim przypadku zakładamy niezależność poszczególnych składowych wektora \(\mathbf{x} = \left(x^1, x^2, ..., x^p \right)^{T}\). W efekcie prowadzi to do następującego wzoru na gęstości prawdopodobieństw \(p(\mathbf{x}|k)\):

\(p(\mathbf{x}|k) = \prod_{i=1}^{p}p(x^i|k) \)

Po zainstalowaniu biblioteki e1071 do zestawu powyższych funkcji dochodzi również naiveBayes() realizująca naiwny klasyfikator Bayesowski. Opcja method="naiveBayes" użyta w funkcji drawparti() wykreśli granice klas dla naiwnego klasyfikatora Bayesowskiego.

library(e1071)

# Klasyfikator naiwnego Bayesa
data.nb <- naiveBayes(class ~ x + y, data)

# Przewidywanie klas za pomocą klasyfikatora naiwnego Bayesa
# domyślnie zwraca klasy
data.nb.pred <- predict(data.nb, data)

# opcja method = "raw" daje prawdopodobieństwa a posteriori
data.nb.pred <- predict(data.nb, data, type = "raw")

# Rysowanie prostej rozdzielającej, punktów etc
with(data, drawparti(class, x, y, method = "naiveBayes", xlab = "X", ylab = "Y", font = 2))

Rysunek 2.8

WIELE KLAS

Nic nie stoi na przeszkodzie, aby powyższe rozważania zastosować do przykładu wielu klas. Poniżej przykładowy kod dla 3 klas

S1 <- matrix(c(4, 0, 0, 4), 2, 2)
S2 <- matrix(c(2, 0, 0, 2), 2, 2)
S3 <- matrix(c(2, 1, 1, 2), 2, 2)

mt1 <- c(-3, -1)
mt2 <- c(2, 2)
mt3 <- c(-2, 3)

n1 <- 40
n2 <- 30
n3 <- 30

X1 <- mvrnorm(n1, mt1, S1)
X2 <- mvrnorm(n2, mt2, S2)
X3 <- mvrnorm(n3, mt3, S3)

X1 <- data.frame(X1); colnames(X1) <- c("x", "y")
X2 <- data.frame(X2); colnames(X2) <- c("x", "y")
X3 <- data.frame(X3); colnames(X3) <- c("x", "y")
X1$class <- 1; X2$class <- 2; X3$class <- 3

plot(X1$x, X1$y, ylim = c(-8,8), xlim = c(-8,8), pch = 19, col = "blue", font = 2)
abline(v = 0, h = 0, col = "gray")
points(X2$x, X2$y, pch = 19, col = "orange")
points(X3$x, X3$y, pch = 19, col = "darkgreen")

Rysunek 2.9

oraz otrzymane rysunki dla metod LDA i QDA

data <- rbind(X1, X2, X3); data$class <- factor(data$class)

# Rysunek dla LDA
with(data, drawparti(class, x, y, xlab = "X", ylab = "Y", font = 2))

# Rysunek dla QDA
with(data, drawparti(class, x, y, method = "qda", xlab = "X", ylab = "Y", font = 2))

Rysunek 2.10

Rysunek 2.11