Zakładamy, że Y = {1, 0}, f 1(~ x) = f (~ x|Y = 1), f 0(~ x) = f (~ x|Y = 0) są gęstościami p-wymiarowego rozkładu normalnego (gaussowskiego)

Klasyfikatory Gaussowskie

Marcin Orchel

1 Wstęp

Wykorzystanie strategii estymacji funkcji gęstości i przyjęcie modelu parametrycznego dla gęstości, tj. znana jest postać gęstości z wyjątkiem tkwiących w niej parametrów.

Zakładamy, że Y = {1, 0}, f ₁ (~ x) = f (~ x|Y = 1), f ₀ (~ x) = f (~ x|Y = 0) są gęstościami p-wymiarowego rozkładu normalnego (gaussowskiego)

f k (~ x) = 1

(2π) ^p/2 |Σ _k | ^1/2 exp



− 1

2 (~ x − ~ µ k ) ^T Σ ⁻¹ _k (~ x − ~ µ k )



(1) Wtedy

X|Y = 1 ∼ N p (µ ₁ , Σ 1 ) (2)

X|Y = 0 ∼ N p (µ ₀ , Σ 0 ) (3)

Twierdzenie 1.1. Jeżeli X|Y = 1 ∼ N _p (µ ₁ , Σ 1 ) i X|Y = 0 ∼ N _p (µ ₀ , Σ 0 ) to klasyfika- tora bayesowski ma postać

d B (~ x) =







1 jeśli r ₁ ² < r ² ₀ + 2 ln ^{ π} _π

0

 + ln ^{ |Σ} _|Σ

^|

1

|



0 poza tym (4)

gdzie

r _i ² = (~ x − ~ µ _i ) ⁰ Σ ⁻¹ _i (~ x − ~ µ _i ) (5) dla i = 0, 1. r _i to odległość Mahalanobisa.

W równoważnej postaci można go zapisać jako d _B (~ x) = arg max

k δ _k (~ x) (6)

gdzie

δ k (~ x) = − 1

2 ln |Σ _k | − 1

2 (~ x − ~ µ k ) ⁰ Σ ⁻¹ _k (~ x − ~ µ k ) + ln π _k (7) Funkcja δ _k nazywa się kwadratową funkcją klasyfikującą grupy k. Procedura klasyfikacji to kwadratowa analiza dyskryminacyjna (QDA).

Nie znamy parametrów π ₁ , π ₀ , µ ₁ , µ ₀ , Σ ₁ , Σ ₀ . Zastępujemy je estymatorami π ˆ 1 = 1

n

n

X

i=1

Y i (8)

π ˆ 0 = 1 n

n

X

i=1

(1 − Y _i ) (9)

µ ˆ ₁ = 1 n ₁

X

Y

X ~ _i (10)

µ ˆ 0 = 1 n 0

X

Y

=0

X ~ i (11)

Σ ˆ ₁ = S ₁ = 1 n ₁ − 1

X

Y

=1

 X ~ _i − ~ µ _i ^{ } X ~ _i − ~ µ ₁ ^{ 0} (12)

Σ ˆ ₀ = S ₀ = 1 n 0 − 1

X

Y

=1

 X ~ _i − ~ µ ₀ ^{ } X ~ _i − ~ µ ₀ ^{ 0} (13)

gdzie

n ₁ =

n

X

i=1

Y _i (14)

n 0 =

n

X

i=1

1 − Y _i (15)

Przy założeniu równości macierzy kowariancji w dwóch grupach, tzn jeśli założymy Σ ₀ = Σ ₁ = Σ klasyfikator bayesowski ma postać

d _B (~ x) = arg max

k δ _k (~ x) (16)

gdzie

δ k (~ x) = ~ x ⁰ Σ ⁻¹ µ ~ k − 1 2

µ ~ ⁰ _k Σ ⁻¹ µ ~ k + ln π _k (17) gdzie k = 1, 0. Funkcja δ _k (~ x) jest funkcją liniową i nosi nazwę liniowej funkcji klasyfiku- jącej grupy k. Powierzchnia rozdzielająca

{~ x : δ ₁ (~ x) = δ ₀ (~ x)} (18) δ ₁₀ (~ x) = {~ x : δ ₁ (~ x) − δ ₀ (~ x)} = 0 (19) jest hiperpłaszczyzną. Jest to liniowa funkcja dyskryminacyjna. Jest ona postaci

σ ₁₀ (~ x) =



~ x − 1

2 (µ ₁ + µ ₀ )

 0

Σ ⁻¹ (µ ₁ − µ ₀ ) + ln

 π ₁ π 0



(20) Obserwację ~ x klasyfikujemy do 1, gdy δ ₁₀ (x) > 0. Procedura klasyfikacji za pomocą liniowej funkcji klasyfikującej nosi nazwę liniowej analizy dyskryminacyjnej (LDA). Nie- znane parametry estymujemy podobnie jak wyżej, z tym, że

Σ = S = ˆ 1

n 1 + n ₀ − 2 ((n ₁ − 1) S ₁ + (n ₁ − 1) S ₀ ) (21)

1.1 Klasyfikator Bayesa dla rozkładów normalnych przy założeniu ta- kich samych macierzy kowariancji

Mamy dane dwa wielowymiarowe rozkłady normalne.

f (~ x|k) = 1

(2π) ^p/2 |Σ| ^1/2 exp



− 1

2 (x − ~ m k ) ^T Σ ⁻¹ (x − ~ m k )



(22) dla k = −1 lub k = 1. Klasyfikator Bayesa maksymalizuje wartość f po obu klasach.

Zobaczmy jak wygląda granica decyzyjna

ln π ₂ f (~ x|Y = 1)

π 1 f (~ x|Y = −1) = 0 (23)

gdzie π ₂ = P (Y = 1), zaś π ₁ = P (Y = −1). A zatem ln π ₂

π 1

− 1

2 ( ~ m 2 − ~ m 1 ) ^T Σ ⁻¹ ( ~ m 1 + ~ m 2 ) + ( ~ m 2 − ~ m 1 ) ^T Σ ⁻¹ ~ x = 0 (24) Widzimy, że jest to funkcja liniowa. Zauważmy, że współczynnik kierunkowy oraz wyraz wolny jest identyczny jak w klasyfikatorze Fishera (przy założeniu, że π ₂ = π ₁ ).

Zwróćmy uwagę na to, że metoda ta estymuje gęstości prawd. P (X|Y = 1), zaś metoda regresji logistycznej estymowała P (Y = 1|X).

1.2 Klasyfikator Bayesa dla rozkładów normalnych przy założeniu róż- nych macierzy kowariancji

Reguła klasyfikacyjna ma postać ln π ₂

π 1

+ 1

2 ln |Σ ₁ |

|Σ ₂ | +x ^{T } Σ ⁻¹ ₂ m ~ ₂ − Σ ⁻¹ ₁ m ~ ₁ ^ − 1

2 x ^{T } Σ ⁻¹ ₂ − Σ ⁻¹ ₁ ^ x− 1

2 m ~ ₂ ^T Σ ⁻¹ ₂ m ~ ₂ + 1

2 m ^T ₁ Σ ⁻¹ m ~ ₁ = 0 (25)

2 Zadania

2.1 Zadania na 3.0

Napisać skrypt w R. W skrypcie

• dla wygenerowanych danych dwuwymiarowych dwóch klas z rozkładów normalnych wykonać klasyfikację lda i qda (macierz kowariancji i wartości średnie obliczone na podstawie danych)

• przetestować klasyfikację dla różnych wartości macierzy kowariancji

• wyświetlić na jednym wykresie punkty, granicę decyzyjną

• wyświetlać na konsoli szybkość wyszukania rozwiązania i jakość mierzoną za po-

mocą średniego błędu klasyfikacji na zbiorze testowym

• wyświetlić punkty, granicę decyzyjną oraz dwie funkcje gęstości prawdopodobień- stwa na wykresie trójwymiarowym

• Dodać komentarz do skryptu opisujący krótko na czym polegają użyte metody oraz wnioski z badań.

Wskazówki do R:

• dane możemy wygenerować za pomocą https://stat.ethz.ch/R-manual/R-devel/

library/MASS/html/mvrnorm.html

• https://stat.ethz.ch/R-manual/R-devel/library/MASS/html/lda.html

• https://stat.ethz.ch/R-manual/R-devel/library/MASS/html/predict.lda.

html

• http://www.statmethods.net/advstats/discriminant.html

• wyraz wolny b nie jest zwracany przez polecenie lda, więc należy go obliczyć, np.

w ten sposób, aby znalezione rozwiązanie przechodziło przez punkt średni.

• polecenie lda działa niepoprawnie dla współliniowych cech

• https://stat.ethz.ch/R-manual/R-devel/library/MASS/html/qda.html

• https://stat.ethz.ch/R-manual/R-devel/library/MASS/html/predict.qda.

html

• http://www.statmethods.net/advstats/discriminant.html

• do rysowania polecenie contour Wskazówki do Matlaba

• https://www.mathworks.com/help/stats/discriminant-analysis.html

• https://www.mathworks.com/help/stats/fitcdiscr.html, deprecated https:

//www.mathworks.com/help/stats/classificationdiscriminant.fit.html, https:

//www.mathworks.com/help/stats/classify.html Przykłady

• uruchomienie lda require("MASS")

sigma <- matrix(c(10,3,3,2),2,2) sigma2 <- matrix(c(10,1,1,5),2,2)

points1<-mvrnorm(n = 1000, c(-1, -1), sigma)

points2<-mvrnorm(n = 1000, c(2, 1), sigma2)

points1b<-cbind(points1, rep(1, 1000))

points2b<-cbind(points2, rep(-1, 1000)) rbind(points1b, points2b)

allPoints<-rbind(points1b, points2b) as.data.frame(allPoints)

allPointsDataFrame<-as.data.frame(allPoints) c <- allPointsDataFrame$V3

w <- allPointsDataFrame$V1 l <- allPointsDataFrame$V2

lda.allPointsDataFrame <- lda(c~w+l, data=allPointsDataFrame)

• uruchomienie lda z wyświetleniem granicy require("MASS")

sigma <- matrix(c(10,3,3,2),2,2) sigma2 <- matrix(c(10,1,1,5),2,2)

points1<-mvrnorm(n = 1000, c(-1, -1), sigma) points2<-mvrnorm(n = 1000, c(2, 1), sigma2) points1b<-cbind(points1, rep(1, 1000)) points2b<-cbind(points2, rep(-1, 1000)) rbind(points1b, points2b)

allPoints<-rbind(points1b, points2b) as.data.frame(allPoints)

allPointsDataFrame<-as.data.frame(allPoints) c <- allPointsDataFrame$V3

w <- allPointsDataFrame$V1 l <- allPointsDataFrame$V2

lda.allPointsDataFrame <- lda(c~w+l, data=allPointsDataFrame) myallPointsDataFrame<- matrix(nrow=2000, ncol=2)

myallPointsDataFrame[,1] <- w myallPointsDataFrame[,2] <- l class<- c

min.w <- min(w) max.w <- max(w) min.l <- min(l) max.l <- max(l)

x<- seq(min.w, max.w, length=100) y<- seq(min.l, max.l, length=100)

myallPointsDataFrameT<-expand.grid(w=x, l=y) n <- length(x)

y.pred <- predict(lda.allPointsDataFrame, as.data.frame(myallPointsDataFrame))

error <- 1 - sum(y.pred$class==class)/length(class)

y.pred <- predict(lda.allPointsDataFrame, myallPointsDataFrameT)

plot(w, l, type="n", xlim=c(min.w, max.w), ylim=c(min.l, max.l))

text(w, l, substr(c, 1, 1))

z <- y.pred$posterior[,1] - y.pred$posterior[,2]

contour(x, y, matrix(z, n), add=TRUE, levels=0.5, labcex = 0, drawlabels = FALSE)

• obliczanie błędu na zbiorze testowym require("MASS")

sigma <- matrix(c(1000,1,3,2),2,2) sigma2 <- matrix(c(1000,1,1,5),2,2)

allPoints1<-mvrnorm(n = 30, c(-2, -1), sigma) allPoints2<-mvrnorm(n = 30, c(2, 1), sigma2) points1b<-cbind(points1, rep(1, 10))

points2b<-cbind(points2, rep(-1, 10)) points1 <- allPoints1[1:10,]

points2 <- allPoints2[1:10,]

points1Testing <- allPoints1[11:30,]

points2Testing <- allPoints2[11:30,]

allTestingPoints<-rbind(points1Testing, points2Testing) points1b<-cbind(points1, rep(1, 10))

points2b<-cbind(points2, rep(-1, 10)) rbind(points1b, points2b)

allPoints<-rbind(points1b, points2b) as.data.frame(allPoints)

allPointsDataFrame<-as.data.frame(allPoints) c <- allPointsDataFrame$V3

classTesting<-append(rep(1, 20), rep(-1, 20)) w <- allPointsDataFrame$V1

l <- allPointsDataFrame$V2

lda.allPointsDataFrame <- lda(c~w+l, data=allPointsDataFrame) allTestingPointsDateFrame <- as.data.frame(allTestingPoints) colnames(allTestingPointsDateFrame) <- c("w", "l")

y.pred <- predict(lda.allPointsDataFrame, allTestingPointsDateFrame) errorTesting <- 1 -

sum(y.pred$class==classTesting)/length(classTesting) myallPointsDataFrame<- matrix(nrow=20, ncol=2)

myallPointsDataFrame[,1] <- w myallPointsDataFrame[,2] <- l class<- c

min.w <- min(w)

max.w <- max(w)

min.l <- min(l)

max.l <- max(l)

x<- seq(min.w, max.w, length=100) y<- seq(min.l, max.l, length=100)

myallPointsDataFrameT<-expand.grid(w=x, l=y) n <- length(x)

y.pred <- predict(lda.allPointsDataFrame,

as.data.frame(myallPointsDataFrame)) error <- 1 - sum(y.pred$class==class)/length(class)

y.pred <- predict(lda.allPointsDataFrame, myallPointsDataFrameT) plot(w, l, type="n", xlim=c(min.w, max.w), ylim=c(min.l, max.l)) text(w, l, substr(c, 1, 1))

z <- y.pred$posterior[,1] - y.pred$posterior[,2]

contour(x, y, matrix(z, n), add=TRUE, levels=0.5, labcex = 0, drawlabels = FALSE)

• rysowanie funkcji 3d różnicy prawdopodobieństw, do poprzedniego przykładu do- dać

persp(x, y, matrix(z, n), theta=150, phi=30, col="lightblue", xlab="x", ylab="y", zlab="z")

• wyświetlanie gęstości prawdopodobieństwa

par(mfrow=c(2,2), pty="m", mar=c(4,2,2,1))

kde2dResult1 <- kde2d(points1[,1], points1[,2], n=25, h=c(3,3), lims=c(-6, 8, -4, 6))

kde2dResult2 <- kde2d(points2[,1], points2[,2], n=25, h=c(3,3), lims=c(-6, 8, -4, 6))

persp(kde2dResult1, phi=20, theta=10, d=5, xlab=c("X1"), ylab=c("X2"))

persp(kde2dResult2, phi=20, theta=10, d=5, xlab=c("X1"), ylab=c("X2"))

• wyświetlanie konturów gęstości prawdopodobieństwa contour(kde2dResult1, nlevels=15, col="red")

contour(kde2dResult2, nlevels=15, col="blue", add=TRUE)

• przykład dla lda, wyświetlanie granicy decyzyjnej names(iris) <- c("d", "sd", "w", "l", "c") c <- iris$c

w <- iris$w l <- iris$l

lda.iris <- lda(c~w+l, data=iris)

myiris<- matrix(nrow=150, ncol=2)

myiris[,1] <- w myiris[,2] <- l class <- c min.w <- min(w) max.w <- max(w) min.l <- min(l) max.l <- max(l)

x<- seq(min.w, max.w, length=100) y<- seq(min.l, max.l, length=100) myirisT<-expand.grid(w=x, l=y) n <- length(x)

y.pred <- predict(lda.iris, as.data.frame(myiris)) error <- 1 - sum(y.pred$class==class)/length(class) y.pred <- predict(lda.iris, myirisT)

plot(w, l, type="n", xlim=c(min.w, max.w), ylim=c(min.l, max.l)) text(w, l, substr(c, 1, 1))

z <- y.pred$posterior[,1] - pmax(y.pred$posterior[,3], y.pred$posterior[,2])

contour(x, y, matrix(z, n), add=TRUE, levels=0.5, labcex = 0, drawlabels = FALSE)

• przykład dla qda, wyświetlanie granicy decyzyjnej names(iris) <- c("d", "sd", "w", "l", "c") c <- iris$c

w <- iris$w l <- iris$l

qda.iris <- qda(c~w+l, data=iris) print(qda.iris)

myiris<- matrix(nrow=150, ncol=2) myiris[,1] <- w

myiris[,2] <- l class <- c min.w <- min(w) max.w <- max(w) min.l <- min(l) max.l <- max(l)

x<- seq(min.w, max.w, length=100) y<- seq(min.l, max.l, length=100) myirisT<-expand.grid(w=x, l=y) n <- length(x)

y.pred <- predict(qda.iris, as.data.frame(myiris))

error <- 1 - sum(y.pred$class==class)/length(class)

print("Error is")

print(error)

y.pred <- predict(qda.iris, myirisT)

plot(w, l, type="n", xlim=c(min.w, max.w), ylim=c(min.l, max.l)) text(w, l, substr(c, 1, 1))

z <- y.pred$posterior[,1] - pmax(y.pred$posterior[,3], y.pred$posterior[,2])

contour(x, y, matrix(z, n), add=TRUE, levels=0.5, labcex = 0, drawlabels = FALSE)

2.2 Zadania na 4.0

• na wykresie trójwymiarowym zaznacz funkcje gęstości dla danych z zadania na 3.0 (macierz kowariancji i wartości średnie obliczone z danych), dane treningowe i klasyfikator qda

• dla wygenerowanych danych trójwymiarowych dwóch klas z rozkładów normal- nych zaznacz na wykresie trójwymiarowym dane treningowe i klasyfikator qda (z macierzą kowariancji i wartościami średnimi obliczonymi z danych)

• oblicz błąd klasyfikacji na zbiorze testowym dla tego klasyfikatora 2.3 Zadania na 5.0

• wykonać klasyfikację qda dla wybranych danych wielowymiarowych ze strony uci,

porównać jakość klasyfikacji z klasyfikatorem Fishera i opartym na regresji logi-

stycznej na danych testowych