W tym raporcie przedstawimy istotność założeń algorytmu k-średnich oraz jedno z zastosowań tego algorytmu

library(car)

library(ClusterR)

#zaladowanie danych z pliku csv

dane <- read.table(file = './dane_zad1.csv',sep=';',dec=',',header = TRUE) 
# przedstawienie danych na wykresie

plot(dane, main='Wizualizacja surowych danych')

Wyraźnie widać zależność obu zmiennych \(V1\) i \(V2\). Rozmieszczenie punktów charakteryzuje się formowaniem dwóch skupisk, jednego centralnego, kołowego skupiska oraz otaczającego go pierścienia punktów.

  1. Poglądowo, tak jak zostało wcześniej wspomniane, można określić, że dane dzielą się na dwie grupy, koło oraz pierścień wokół niego.

  2. Teraz określimy optymalną liczbę klastrów, korzystając z wykresu osypiska, gdzie na osi \(y\) zaznaczymy całkowitą wewnątrzgrupową sumę kwadratów odległości (total within-clusters sum of squares), a na osi \(x\) mamy liczbę klastrów. W tym celu w pętli wykonujemy algorytm k-średnich dla różnej liczby klastrów. Otrzymane rezultaty zaznaczamy na wykresie oraz porównujemy wykres z wykresem wygenerowanym za pomocą funkcji

# przeprowadzenie algorytmu kmeans
# a) pogladowo 2 klastry

# dobranie odpowiedniej liczby klastrow
# b) metoda lokcia
k <- 10 # Maksymalna liczba klastrów
t.wss<-c()
#wss <- (nrow(dane)-1)*sum(apply(dane,2,var))
for (i in 1:k) t.wss[i] <- kmeans(dane, centers=i)$tot.withinss

plot(1:k, t.wss,
     type="b", pch = 16, 
     xlab="liczba klastrów",
     ylab="Łączna suma wariancji klastrów",
     main = 'Wykres osypiska',
     col="blue"
)
abline(v = 6, lty =2,col='red')
grid()

# porownanie z funkcja Optimal_Cluster_KMeans

Optimal_Clusters_KMeans(dane,max_clusters = 10,criterion = 'WCSSE')

##  [1] 7737.2112 5522.7577 3688.1585 2391.7826 1731.8769 1280.7342 1025.3866
##  [8]  881.1016  776.3808  618.1568

Oba wykresy są podobne, różnią się nieco wartościami dla poszczególnych klastrów, co jest spowodowane faktem, że algorytm za każdym razem losowo wybiera początkowe centroidy, co wpływa na ostateczne wyniki. Jednak z wykresów można wyciągnąć te same wnioski, że odpowiednią liczbą klastrów jest 6, potem tot.withinss nie maleje już drastycznie.

Teraz wykonamy grupowanie metodą k-średnich, najpierw dzieląc na 2, potem na 6 klastrów. W R służy do tego funkcja kmeans(), ustawiamy liczbę klastrów, zostawiamy domyślny algorytm “Hartigana-Wonga”. Rezultaty naszego grupowanie przedstawiamy na wykresie, kwadracikiem zaznaczamy centroidy klastrów.

# algorytm k srednich dla 2 klastrow
podzial_a <- kmeans(dane, centers = 2)
wektor_podzialu_a <- as.factor(podzial_a$cluster) #wektor podzialu

#wizualizacja podzialu
plot(dane[,1], dane[,2], col =wektor_podzialu_a, pch = 19,xlab = 'v1',
     ylab = 'v2',main = 'Podzial na 2 klastry')
# naniesienie srodkow
points(podzial_a$centers[,1],podzial_a$centers[,2],col=c('grey','red'),pch=46,cex=13)

# algorytm k srednich dla 6 klastrow
podzial_b <- kmeans(dane, centers = 6)
wektor_podzialu_b <- as.factor(podzial_b$cluster)

#wizualizacja podzialu
plot(dane[,1], dane[,2], col =wektor_podzialu_b, pch = 19,xlab = 'v1',
     ylab = 'v2',main = 'Podzial na 6 klastrow')
# naniesienie srodkow
points(podzial_b$centers[,1],podzial_b$centers[,2],
       col=c('grey','red','green','blue','steelblue','pink'),pch=46,cex=13)

W obu przypadkach widzimy, że grupowanie nie zostało przeprowadzony tak jak byśmy chcieli, grupy nie odzwierciedlają naturalnie formujących się w danych skupisk. To ma związek z tym, że algorytmu korzysta z odległości euklidesowej, przez co buduje sferyczne skupiska. W naszych danych mamy koło i pierścień wokół niego, zatem nie jesteśmy w stanie za pomocą sferycznych skupisk rozdzielić tych dwóch grup.

Rozwiązaniem naszego problemu może być zamiana zmiennych na współrzędne biegunowe. W tym celu, korzystając ze wzorów \(r=\sqrt{v1^2+v2^2}\) i \(\theta=\arctan{(\frac{v2}{v1})}\), przekształcamy nasze zmienny. Dane po przekształceniu wyglądają tak:

# zmiana wspolrzednych na wspolrzedne biegunowe
r <- sqrt(dane[,1]^2+dane[,2]^2)
theta <- atan(dane[,2]/dane[,1])

#tabela z nowymi wspolrzednymi
dane2 <- cbind(r,theta)
plot(dane2,main = "Rozproszenie danych w nowym ukladzie")

Możemy zauważyć 2 osobne zbiory, które to będzie można zgrupować korzystając z algorytmu k-średnich. Wykonujemy wykres osypiska, z którego można odczytać, że optymalną liczbą klastrów są dwa klastry.

# znalezienie optymalnej liczby klastrow
Optimal_Clusters_KMeans(dane2,max_clusters = 10,criterion = 'WCSSE')

##  [1] 2700.22160  462.16997  335.90231  210.68596  184.12120  144.73596
##  [7]  134.04970  113.12951   94.62957   83.78474
# algorytm k srednich dla 2 klastrow
podzial <- kmeans(dane2, centers = 2)
wektor_podzialu <- as.factor(podzial$cluster)

#wspolrzedne srodka
x_c<-podzial$centers[,1]*cos(podzial$centers[,2])
y_c<-podzial$centers[,1]*sin(podzial$centers[,2])

#wizualizacja podzialu
plot(dane[,1], dane[,2], col =wektor_podzialu, pch = 19,xlab = 'v1',
     ylab = 'v2',main = 'Podzial na 2 klastry')
points(x_c,y_c,col=c('grey','red'),pch=46,cex=13)

Następnie przeprowadzamy algorytm k-średnich i wydobywamy etykiety przyporządkowujące punkty do grup oraz przedstawiamy rezultaty na wyjściowych danych.

Można zauważyć, że tym razem algorytm pogrupował dane zgodnie z naszymi oczekiwaniami.

Klasteryzacja hierarchiczna

Klasteryzacja metodą pojedynczego wiązania i metodą pełnego wiązania należą do grupy algorytmów hierarchicznych aglomeracyjnych. W obu przypadkach, na początku wszystkie punkty traktujemy jako osobne klastry, w kolejnych iteracjach łączymy klastry na zasadzie podobieństwa (w sensie odległości), aż do momentu, gdy wszystkie punkty będą należeć do jednego klastra. Różnica między tymi metodami polega na sposobie mierzenia odległości między klastrami. W algorytmie pojedynczego wiązania, odległość między dwoma klastrami jest minimalną odległością między elementami klastrów. Natomiast w algorytmie całkowitego wiązania, jako odległość bierzemy maksymalną odległość między elementami badanych klastrów.

Na początku tworzymy macierz odległości między wszystkimi punktami (w R służy do tego funkcja dist()}. Następnie szukamy najbardziej podobnych klastrów i grupujemy je, wyliczamy odległości i aktualizujemy naszą macierz. Powtarzamy to iteracyjnie, aż do momentu, gdy wszystkie obiekty znajdujące się w jednym klastrze.

Wykonujemy algorytm w R, korzystając z funkcji hclust(), ustawiając metodę najpierw jako “single” dla pojedynczego wiązania, następnie “complete” dla pełnego wiązania. Rysujemy dendrogramy.

#klasteryzacja hierarchiczna dla pojedynczego wiazania
#funkcja dist() tworzy macierz odleglosci
clust_s<-hclust(dist(dane), method = 'single')
plot(clust_s) #dendogram

label1 <-cutree(clust_s,k=2) #podzial na dwa klastry, tniemy dendogram


#klasteryzacja hierarchiczna dla calkowitego wiazania
clust_c<-hclust(dist(dane), method = 'complete')
plot(clust_c) #dendogram

label2 <-cutree(clust_c,k=2) #podzial na dwa klastry, tniemy dendogram

Funkcją cutree() odcinamy dendrogramy dla ustalonej liczby klastrów, u nas dla \(k=2\) i pobieramy wektor etykiet. Przedstawiamy rezultaty oby grupowań.

plot(dane[,1], dane[,2], col =label1, pch = 19,xlab = 'v1',
     ylab = 'v2',main = 'Podzial na 2 klastry algorytmem pojedynczego wiazania')

plot(dane[,1], dane[,2], col =label2, pch = 19,xlab = 'v1',
     ylab = 'v2',main = 'Podzial na 2 klastry algorytmem calkowitego wiazania')

Widzimy, że algorytm pojedyńczego wiązania poprawnie podzielił dane. Metoda pełnego wiązania nie poradziła sobie w tym przypadku, ta metoda tworzy bardziej zwarte skupiska, przez co w naszym wypadku nie zadziałała, tak jak byśmy chcieli.

PCA i klasteryzacja

Wczytujemy dane z pakietu DAAG, do analizy wyciągamy odpowiednie zmienne oraz zapisujemy wektor przechowujący płeć obiektów.

#install.packages('DAAG')
library(DAAG)
head(ais)
##    rcc wcc   hc   hg ferr   bmi   ssf pcBfat   lbm    ht   wt sex  sport
## 1 3.96 7.5 37.5 12.3   60 20.56 109.1  19.75 63.32 195.9 78.9   f B_Ball
## 2 4.41 8.3 38.2 12.7   68 20.67 102.8  21.30 58.55 189.7 74.4   f B_Ball
## 3 4.14 5.0 36.4 11.6   21 21.86 104.6  19.88 55.36 177.8 69.1   f B_Ball
## 4 4.11 5.3 37.3 12.6   69 21.88 126.4  23.66 57.18 185.0 74.9   f B_Ball
## 5 4.45 6.8 41.5 14.0   29 18.96  80.3  17.64 53.20 184.6 64.6   f B_Ball
## 6 4.10 4.4 37.4 12.5   42 21.04  75.2  15.58 53.77 174.0 63.7   f B_Ball
dane <- ais[,1:11] #wybranie danych
plec <- ais[,12] #wektor etykiet

Korzystając z funkcji prcomp() wykonujemy analizę PCA na nieprzeskalowanych danych i przedstawiamy rzutowanie na pierwsze trzy składowe z rozróżnieniem na płeć.

library(factoextra)
library(stats)
library(ggpubr)
library(plotly)
#implementacja PCA przy wykorzystaniu funkcji prcomp(), dane nieprzeskalowane
dane.pca <- prcomp(dane,center = F,scale = F,retx = T)

# wykres rozproszenia, projekcja PC1 i PC2
wyk_diag12<-fviz_pca_ind(dane.pca,
                         axes= c(1,2),
                         col.ind = plec, 
                         habillage = plec,
                         geom = c("point"),
                         palette = c(1,  2),
                         legend.title = "Sex",
                         repel = T
)
ggpubr::ggpar(wyk_diag12, main = 'Projekcja obserwacji na PC1 i PC2.',
              xlab = "PC1", ylab = "PC2")

# wykres rozproszenia, projekcja PC1 i PC3
wyk_diag13<-fviz_pca_ind(dane.pca,
                         axes= c(1,3),
                         col.ind = plec, 
                         habillage = plec,
                         geom = c("point"),
                         palette = c(1, 2),
                         legend.title = "Sex",
                         repel = T
)
ggpubr::ggpar(wyk_diag13, main = 'Projekcja obserwacji na PC1 i PC3.',
              xlab = "PC1", ylab = "PC3")

# wykres rozproszenia PC2 i PC3
wyk_diag23<-fviz_pca_ind(dane.pca,
                         axes= c(2,3),
                         col.ind = plec, 
                         habillage = plec,
                         geom = c("point"),
                         palette = c(1,  2),
                         legend.title = "Sex",
                         repel = T
)
ggpubr::ggpar(wyk_diag23, main = 'Projekcja obserwacji na PC1 i PC3.',
              xlab = "PC2", ylab = "PC3")

# wykres 3d rozproszenia, projekcja PC1, PC2 i PC3

scatter3d <- plot_ly(x = dane.pca$x[,1], y = dane.pca$x[,2], z = dane.pca$x[,3],
                     color = plec, colors = c("#00AFBB",  "#FC4E07")) %>%
  
  add_markers(size = 12)

tytul <- "Projekcja na PC1,PC2 i PC3"
scatter3d <- scatter3d %>%
  
  layout(
    
    title = tytul,
    scene = list(bgcolor = "#ffffff")
    
  )
scatter3d

Na rzutowaniu PC1, PC2 trudno dostrzec jakieś grupowanie, grupy są ze sobą raczej wymieszane. Na pozostałych projekcjach grupy formują się już wyraźniej, chociaż niekoniecznie w sferyczne grupy.

Wykonujemy algorytm k-średnich dla pierwszych dwóch składowych, dzielimy dane na dwa klastry, ponieważ chcemy dokonać dyskryminacji ze względu na płeć. Kolorem (czerwony i czarny) rozróżniamy grupy, które powstały za pomocą algorytmu k-średnich, natomiast kształtem (kółko i trójkąt) zaznaczamy rzeczywistą płeć. Legenda

#klasteryzacja pierwszych dwoch skladowych
set.seed(6541)
kM <- kmeans(x=cbind(dane.pca$x[,1],dane.pca$x[,2]), centers = 2)
wek_kM <- as.factor(kM$cluster)

#plec 1-kolko-kobieta, 2-trojkat-mezczyzna
#wek_kM 1-czarny, 2- czerwony
# czarne kolko - Tf, czerwone kolko -Ff, czerwony trojkat -Tm
plot(dane.pca$x[,1:2],col=wek_kM,pch=as.numeric(plec),xlab = 'PC1',ylab = 'PC2',
            main ='Rezultat przyporzadkowania vs. rzeczywistosc')
legend(-200,130,legend=c('Tf','Ff','Tm','Fm'),col=c(1,2,2,1),pch=c(1,1,2,2))

Widzimy że jest wiele obserwacji, które zostały nieprawidłowo rozszyfrowane. Algorytm szczególnie wielu mężczyzn nieprawidłowo zgrupował jako kobiety.

Następnie wykonujemy analizę PCA jeszcze raz, tym razem na przeskalowanych danych, a na wykresie przedstawiamy rzutowanie pierwszych trzech składowych.

#implementacja PCA przy wykorzystaniu funkcji prcomp(), dane przeskalowane
dane.pca2 <- prcomp(dane,center = T,scale = T,retx = T)

# wykres rozproszenia, projekcja PC1 i PC2
wyk_diag21<-fviz_pca_ind(dane.pca2,
                         axes= c(1,2),
                         col.ind = plec, 
                         habillage = plec,
                         geom = c("point"),
                         palette = c(1,  2),
                         legend.title = "Sex",
                         repel = T
)
ggpubr::ggpar(wyk_diag21, main = 'Projekcja obserwacji na PC1 i PC2.',
              xlab = "PC1", ylab = "PC2")

# wykres rozproszenia, projekcja PC1 i PC3
wyk_diag22<-fviz_pca_ind(dane.pca2,
                         axes= c(1,3),
                         col.ind = plec, 
                         habillage = plec,
                         geom = c("point"),
                         palette = c(1,  2),
                         legend.title = "Sex",
                         repel = T
)
ggpubr::ggpar(wyk_diag22, main = 'Projekcja obserwacji na PC1 i PC3.',
              xlab = "PC1", ylab = "PC3")

# wykres rozproszenia PC2 i PC3
wyk_diag23<-fviz_pca_ind(dane.pca2,
                         axes= c(2,3),
                         col.ind = plec, 
                         habillage = plec,
                         geom = c("point"),
                         palette = c(1,  2),
                         legend.title = "Sex",
                         repel = T
)
ggpubr::ggpar(wyk_diag23, main = 'Projekcja obserwacji na PC2 i PC3.',
              xlab = "PC2", ylab = "PC3")

# wykres 3d rozproszenia, projekcja PC1, PC2 i PC3

scatter3d2 <- plot_ly(x = dane.pca2$x[,1], y = dane.pca2$x[,2], z = dane.pca2$x[,3],
                     color = plec, colors = c("#00AFBB",  "#FC4E07")) %>%
  
  add_markers(size = 12)

tytul <- "Projekcja na PC1,PC2 i PC3"
scatter3d2 <- scatter3d2 %>%
  
  layout(
    
    title = tytul,
    scene = list(bgcolor = "#ffffff")
    
  )
scatter3d2

W tym wypadku na rzutowaniu PC1, PC2 (jak i PC1, PC3) widać dwie odseparowane grupy, dla których można zastosować algorytm k-średnich. Rzut na drugą i trzecia składową nie przedstawia dwóch grup, natomiast do dyskryminacji ze względu na płeć wystarczą dwie pierwsze składowe.

Na wykresie trójwymiarowym również można dostrzec formowanie się grup.

Analogicznie jak wcześniej wykonujemy algorytm k-średnich i zaznaczamy rezultat na wykresie, wraz z rzeczywistymi etykietami płci. Legenda i opis jak wyżej.

#klasteryzacja pierwszych dwoch skladowych
set.seed(6541)
kM2 <- kmeans(x=cbind(dane.pca2$x[,1],dane.pca2$x[,2]), centers = 2)
wek_kM2 <- as.factor(kM2$cluster)

#plec 1-kolko-kobieta, 2-trojkat-mezczyzna
#wek_kM 1-czarny, 2- czerwony
# czarne kolko - Tf, czerwone kolko -Ff, czerwony trojkat -Tm
plot(dane.pca2$x[,1:2],col=wek_kM2,pch=as.numeric(plec),xlim=c(-7,6),xlab = 'PC1',ylab = 'PC2',
     main ='Rezultat przyporzadkowania vs. rzeczywistosc')
legend(-7,5,legend=c('Tf','Ff','Tm','Fm'),col=c(1,2,2,1),pch=c(1,1,2,2))

W tym wypadku widzimy, że przyporządkowanie algorytmu jest dużo lepsze, nadal są obserwacje nieprawidłowo przyporządkowane do płci, ale jest ich znacznie mniej.

Zatem możemy wnioskować, że należy przeprowadzać analizę na przeskalowanych danych.

Zastosowanie k-średnich w kompresji

Zajmiemy się teraz kompresją obrazu przy pomocy algorytmu k-średnich. Wczytujemy obrazek do R, zapisujemy na komputerze oraz odczytujemy kanały odpowiednich kolorów. Najpierw zajmiemy się kompresją obrazka z rozszerzeniem .png.

#install.packages('imager')
library(imager)
#install.packages('OpenImageR')
library(OpenImageR)

#algorytm dla png
fpath <- system.file('extdata/parrots.png',package='imager')
papugi <- readImage(fpath) #odczytanie obrazka
#imageShow(fpath)

#writeImage(papugi,'papugi.png')
#odczytanie kanalow
r <- papugi[,,1]
g <- papugi[,,2]
b <- papugi[,,3]

Papugi zdjecie oryginalne

Najpierw tworzymy siatkę 512x768, która odpowiada za rozmieszczenie pikseli i następnie tworzymy ramkę danych zawierającą te współrzędne pikseli oraz kolory RGB.

#utworzenie siatki ze wspolrzednymi, siatka 768 szer, 512 wys
x<-seq(1,768,1)
y<-seq(512,1,-1)
wymiar<-expand.grid(y=y,x=x)

#zapisanie wspolrzedne + RGB w tabeli
imRGB <- transform(wymiar, R = r[as.matrix(wymiar)], G=g[as.matrix(wymiar)], B=b[as.matrix(wymiar)])
head(imRGB)
##     y x         R         G         B
## 1 512 1 0.0000000 0.0000000 0.0000000
## 2 511 1 0.2392157 0.3333333 0.1137255
## 3 510 1 0.2392157 0.3333333 0.1254902
## 4 509 1 0.2392157 0.3490196 0.1333333
## 5 508 1 0.2627451 0.3686275 0.1568627
## 6 507 1 0.2588235 0.3803922 0.1529412

Dokonujemy algorytmu k-średnich dla \(k=\{2,16,64,256\}\).

#algorytm k srednich dla kanalow dla k={2,16,64,256}
# podmieniam k i puszczam od nowa, rezultaty zapisuje w raporcie
k<-2
alg<-kmeans(imRGB[,c('R','G','B')],centers = k,iter.max = 300)
wek<-alg$cluster

Wektor etykiet zapisujemy do tabeli klastry, a centroidy wraz z numerem klastra zapisujemy w tabeli o nazwie centra.

klastry<-data.frame(nr_klastra=wek)
head(klastry)
##   nr_klastra
## 1          1
## 2          1
## 3          1
## 4          1
## 5          1
## 6          1
#tabela centroid + kod klastra
centra <- data.frame(nr_klastra=(1:k),alg$centers)
head(centra)
##   nr_klastra         R         G         B
## 1          1 0.3884209 0.3504322 0.2324227
## 2          2 0.8003630 0.7190386 0.5334574

Kolejnym krokiem jest zakodowanie, czyli dla każdego piksela przyporządkowujemy numer klastra, do którego należy. Tak utworzoną ramkę zapisujemy jako img_kod.

#tabela kodujaca pixele wedlug nr klastra
img_kod<-data.frame(imRGB[,1:2],klastry)
head(img_kod)
##     y x nr_klastra
## 1 512 1          1
## 2 511 1          1
## 3 510 1          1
## 4 509 1          1
## 5 508 1          1
## 6 507 1          1

Aby zdekodować informacje, musimy połączyć tabele, tak aby nowa zawierała wszystkie etykiety klastrów wraz z przyporządkowanymi im kolorami RGB.

library(tidyverse)
#dekodowanie obrazka, laczenie tabeli z centrami i tabeli ze wspolrzednymi
kod_rgb<-full_join(img_kod,centra,by = "nr_klastra")
kod_rgb<-kod_rgb[,3:6]
head(kod_rgb)
##   nr_klastra         R         G         B
## 1          1 0.3884209 0.3504322 0.2324227
## 2          1 0.3884209 0.3504322 0.2324227
## 3          1 0.3884209 0.3504322 0.2324227
## 4          1 0.3884209 0.3504322 0.2324227
## 5          1 0.3884209 0.3504322 0.2324227
## 6          1 0.3884209 0.3504322 0.2324227

Następnie tworzymy trójwymiarową macierz o wymiarach 512x768x3 i taką tablicę możemy zapisać na dysku funkcją writeImage(), a także otworzyć w programie funkcją OpenImageR().

#utworzenie macierzy z kolorami do odczytu obrazka
macierz<-as.matrix(kod_rgb[,2:4],nrow = 512,ncol=768,byrow = T)
papugi_komp<-array(macierz, dim = c(512,768,3))

#writeImage(flipImage(papugi_komp,mode = 'vertical') ,'papugi_2.png')
#imageShow('./papugi_2.png')

Poniżej przedstawiamy rezultaty kompresji dla \(k=\{2,16,64,256\}\) kolejno.

papuga 2 kolory papuga 16 kolory papuga 64 kolory papuga 256 kolory

Oryginalny obrazek miał 549 KB, a kolejne kompresje mają 13.3 KB, 78.1 KB, 171 KB, 284 KB. Skompresowane obrazki rzeczywiście mają mniejsze rozmiary, jednak oczywiście tracimy pewien poziom szczegółowości - zmniejszamy paletę użytych kolorów.