library(corrplot)

# Baza USArrests sadrzi podatke o hapsenje za nekoliko kriminalnih prekrsaja,
# kao i procenat urbane populacije, za svaku od 50 americkih saveznih drzava
data("USArrests")
data = USArrests
corrplot(cor(USArrests)) #podaci su logicno, jako korelisani

# Nad podacima cemo da uradimo kmeans klasterizaciju, koje podatke nenadgledano,
# odnosno bez ciljne promenljive, deli na k grupa (klastera)
scaled_data <- data.frame(scale(USArrests))
kmeans_model <- kmeans(scaled_data, centers = 5, nstart = 10) #Ponavljamo algoritam 10 puta, bira se samo najbolji ishod
kmeans_model
# Model vraca razne podatke o klasterizaciji: dodeljeni klasteri, sredine
# klastera, unutarklasterna i medjuklasterno rastojanje, itd.

set.seed(216)
wss <- c()
for(k in 1:20)
  wss[k] = kmeans(scaled_data, k, nstart = 10)$tot.withinss #variranjem vrednost nstart dobijamo manje ili vise pravilnu krivu

plot(1:20, wss, type = "o", main = "Unutarklasterna suma kvadrata u zavisnosti od broja klastera") #lakat na k=4

final_model <- kmeans(scaled_data, centers = 4, nstart = 10)
final_model$cluster #dodeljeni klasteri svih americkih drzava

# Mozemo da vidimo, na primer, kako klasteri izgledaju na grafiku 2 atributa
cluster <- final_model$cluster
plot(scaled_data$Murder, scaled_data$UrbanPop, pch = 16, cex = 1.5, col = cluster)
text(scaled_data$Murder, scaled_data$UrbanPop, labels = rownames(scaled_data), pos = 3, cex = 0.8)

# Sada, kada bismo sa ovim podacima radili neku unarsknu validaciju, ili na bilo
# koji nacin delili skup, mogli da uzorkujemo iz klastera nezavisno, kako bismo
# imali jednaku reprezentaciju drzava iz svakog klastera

# Mane kmeans-a:

library(MASS)

# Generisemo podatke iz visedimenzionih normalnih raspodela
set.seed(216)

mean1 <- c(0, 0)
cov_matrix1 <- matrix(c(25,0,0, 1), nrow = 2)
normal1 <- mvrnorm(100, mean1, cov_matrix1)

mean2 <- c(0, 5)
cov_matrix2 <- matrix(c(25, 0, 0, 1), nrow = 2)
normal2 <- mvrnorm(100, mean2, cov_matrix2)

gen_data <- rbind(normal1, normal2)
group <- c(rep('blue',100),rep('red',100))

plot(gen_data[,1], gen_data[,2], pch = 16, asp = 1) #jasno izdvojene elipse
#plot(gen_data[,1], gen_data[,2], col = group, pch = 16, asp = 1) #da vidimo iz kojih raspodela poticu

# Kako radi kmeans?
clusters = kmeans(gen_data, centers = 2, nstart = 10)$cluster
plot(gen_data[,1], gen_data[,2], col = clusters, pch = 16, asp = 1)

# Mozemo, izmedju ostalog, i da pogledamo kakvo k bi kmeans preporucio
wss <- c()
for(k in 1:20)
  wss[k] = kmeans(gen_data, k, nstart = 10)$tot.withinss
plot(1:20, wss, type = "o", main = "Unutarklasterna suma kvadrata u zavisnosti od broja klastera")
# Ne postoji jasan lakat

## Mesavina normalnih raspodela/ EM algoritam

library(mclust)

gmm_model <- Mclust(gen_data, 2)
clusters <- gmm_model$classification
plot(gmm_model, what = 'classification', asp = 1) #gmm ima lep ugradjeni plot metod
gmm_model$parameters #parametri su dosta blizu pravih

# Mesavina normalnih raspodela (Gaussian Mixture Model - GMM) predvidja ne samo
# grupe, vec i verovatnoce pripadanja bilo kojim od grupa
gmm_model$uncertainty #nesigurnost modela u svoju odluku

# Kod koji pravi neprekidnu skalu boja
color_scale = colorRampPalette(c("blue","red")) 
col = color_scale(100)[cut(gmm_model$uncertainty, breaks = 100)]

# Crtamo nesigurnost modela, crvenije tacke su manje sigurne
plot(gen_data[,1], gen_data[,2], col = col, pch = 16, asp = 1)
# Vidimo da je pri ovim podacima model bio jako siguran u odlukama za skoro
# sve tacke

# U R implementaciji imamo mogucnost pravljenja modela za razlicite vrednosti
# broja grupa istovremeno, odnosno moze da se koristi i da automatski
# izabere broj grupa.
# Podrazumevano proverava vrednosti od 1 do 9
Mclust(gen_data)$G #najbolje su 2 grupe, kao sto smo i ocekivali  

# Na bazi USArrests
gmm_model <- Mclust(scaled_data)
gmm_model$G #preporucene su 3 grupe
gmm_model$uncertainty #mali brojevi, prilicno je siguran u svoje odluke
cluster = gmm_model$classification
plot(scaled_data$Murder, scaled_data$UrbanPop, pch = 16, cex = 1.5, col = cluster)
text(scaled_data$Murder, scaled_data$UrbanPop, labels = rownames(scaled_data), pos = 3, cex = 0.8)
# Prakticno su 2 klastera koja smo imali kod kmeans-a stopljena
plot(gmm_model, what = 'classification') #grafici svih kombinacija atributa


## DBSCAN

# Za razliku od prethodna dva modela/algoritma, dbscan ne pokusava da uoci
# centralizovane grupe podataka, vec guste regione proizvoljnog oblika
# Algoritam krece od neke glavne tacke, koju definisemo kao tacku sa barem
# z tacaka u okolini. U njen klaster dodaje sve tacke iz okoline,
# i onda se tako klaster 'siri'.
# U ovom algoritmu ne definisemo broj klastera, vec sirinu okoline - epsilon
# i minimalan broj okolnih komsija za glavne tacke - z

library(factoextra)
library(dbscan)

data(multishapes)
df <- multishapes[,1:2] #treca promenljiva je kategoricka i oznacava grupu
plot(df, asp = 1)
# Ne ocekujemo ni od kmeans-a ni GMM-a da moze da razdvoji dve kruznice sa slike


# Probacemo sa z = 2, generalno je preporuceno da z bude istog reda velicine kao
# broj atributa
z = 2
# Crtamo grafik sortiranih distanci drugog najblizeg komsije za sve tacke u df
kNNdistplot(df, k = z)
abline(h = 0.13)
# 0.13 je otprilike na laktu, pa to biramo za epsilon
eps = 0.13
dbscan_model <- dbscan(df, eps = eps, minPts = z)
clusters <- dbscan_model$cluster
clusters #nule su tacke suma, nisu dodeljene nijednom klasteru
plot(df[,1], df[,2], col = clusters + 1, pch = 16, asp = 1) #clusters + 1, jer inace ne crta tacke suma
# Tacke suma su crne
# DBSCAN je dobro uhvatio pravilnosti u podacima, ali vidimo dosta nepotrebnih
# malih klastera od po 2 tacke
# Ovo mozemo da izbegnemo sa vecom vrednosti z
clusters <- dbscan(df, eps = eps, minPts = 4)$cluster
plot(df[,1], df[,2], col = clusters + 1, pch = 16, asp = 1)
# Sada deluje idealno, ali u prakticnim situacijama, kada imamo visedimenzione
# podatke, ne mozemo da vidimo ovo sa slike
# Kada bismo za z = 4 odredjivali epsilon pravilom lakta:
kNNdistplot(df, k = 4)
abline(h = 0.17) #izgleda dobro
clusters <- dbscan(df, eps = 0.17, minPts = 4)$cluster
plot(df[,1], df[,2], col = clusters + 1, pch = 16, asp = 1)
# Sada su dve kruznice stopljene u isti klaster
# Ako probamo malo manje epsilon, recimo 0.12:
clusters <- dbscan(df, eps = 0.12, minPts = 4)$cluster
plot(df[,1], df[,2], col = clusters + 1, pch = 16, asp = 1)
# Sada imamo dodatan klaster na kruznici
# Ovaj algoritam je generalno jako osetljiv na izbor epsilona, koji nije nimalo
# ocigledan

# Probajmo, primera radi, kmeans nad podacima
wss <- c()
for(k in 1:20)
  wss[k] = kmeans(df, k, nstart = 10)$tot.withinss 
plot(1:20, wss, type = "o", main = "Unutarklasterna suma kvadrata u zavisnosti od broja klastera")
# Deluje da je lakat na k=5, mozda i na 2
k = 5
cluster_kmeans = kmeans(df, centers = k, nstart = 10)$cluster
plot(df[,1], df[,2], col = cluster_kmeans, pch = 16, asp = 1)
k = 2
cluster_kmeans = kmeans(df, centers = k, nstart = 10)$cluster
plot(df[,1], df[,2], col = cluster_kmeans, pch = 16, asp = 1)