library(mlbench)
library(FNN)
library(caret)
library(dplyr)
library(ggplot2)

data(PimaIndiansDiabetes)
data = PimaIndiansDiabetes
N = nrow(data)

# Podsetimo se kako knn model izgleda za razlicite vrednosti k
# Kod sa treceg casa:
X <- data[,c(2,6)]
y <- as.numeric(data$diabetes)
glucose_dense = seq(-5, 205, length.out = 100)
mass_dense = seq(-5, 70, length.out = 100)
df_dense = expand.grid(glucose = glucose_dense, mass = mass_dense)

dense_predictions = knn(train = X, test = df_dense, cl = y, k = 51) #ovde varirati vrednost k
df = df_dense
df$prediction = dense_predictions
ggplot(data = df, aes(x = glucose, y = mass, color = as.factor(prediction))) +
  geom_point()
# Sto je vece k, to je model jednostavniji i manje se preprilagodjava podacima
# Ali sa druge strane jednostavniji model teze uvidja kompleksne obrasce u podacima

# Hocemo da vidimo kako se kompleksnost modela odrazava na njegovo performanse
# na trening i test skupu. Ocekujemo da kompleksniji modeli imaju bolje
# performanse na trening skupu. Sa druge strane, kod jednostavnijih modela je
# manja razlika izmedju performansi na trening i test skupu

# Evaluiramo tacnost knn-a za razlicite vrednosti parametra k (1-100) na 
# trening i test skupu koristeci unakrsnu validaciju
# Izvrsavanje koda treba da traje neko vreme (meni je trebalo 18s), zato sto
# pravimo 100*10 modela
# Za brze izvrsavanje smanjiti broj foldova
ks = 1:100
# Pravimo foldove za unakrsnu validaciju
n_folds = 10
fold_idx = createFolds(data$diabetes, k = n_folds)
folds = lapply(fold_idx, function(idx) data[idx,])

cv_performance_by_k = sapply(ks, function(k) {
  train_accs = c()
  test_accs = c()
  for(i in 1:n_folds)
  {
    train_set = bind_rows(folds[-i])
    test_set = folds[[i]]
    
    #podela na trening i test
    X_train = train_set[,-9]
    y_train = train_set[,9]
    X_test = test_set[,-9]
    y_test = test_set[,9]
    
    #skaliranje
    X_train_means = colMeans(X_train)
    X_train_sds = apply(X_train, 2, sd)
    X_train_scaled = data.frame(scale(X_train, center = X_train_means, scale = X_train_sds))
    X_test_scaled = data.frame(scale(X_test, center = X_train_means, scale = X_train_sds))
    
    #predvidjanje
    y_pred_train = knn(X_train_scaled, X_train_scaled, y_train, k = k)
    y_pred_test = knn(X_train_scaled, X_test_scaled, y_train, k = k)
    
    #tacnost predvidjanja
    train_accs[i] = mean(y_pred_train == y_train)
    test_accs[i] = mean(y_pred_test == y_test)
  }
  return (c(mean(train_accs), mean(test_accs)))
})
plot(ks, cv_performance_by_k[1,], type = "l", col = "red", xlab = "K", ylab = "Accuracy", ylim = c(0.65, 1))
title(main = "Train and test peformance of KNN by K")
lines(ks, cv_performance_by_k[2,], col = "blue")
legend("topright", legend = c("Train Set", "Test Set"), 
       col = c("red", "blue"), lty = 1)

# Kao što smo napomenuli, generalni trend je da performanse opadaju na trening
# skupu, kako se k povecava. Uz to, performanse su uvek bolje na trening skupu.
# Za jako male vrednosti k imamo jasno preprilagodjavanje - performanse su mnogo
# bolje na trening skupu. Za selekciju najbolje vrednosti k gledamo perforanse
# na trening skupu: Vidimo da je idealna neka vrednost k izmedju 20 i 80,
# sve manje ili vece daje jasno losije performanse.
# Formalno, biramo k koje je bilo najbolje na test skupu:
best_k = ks[which.max(cv_performance_by_k[2,])]
best_k
max(cv_performance_by_k[2,])
abline(v = best_k)
abline(h = max(cv_performance_by_k[2,]))
# U ovom slucaju smo odvojili test skup kako bismo odabrali najbolju vrednost
# hiperparametre (ili najbolji model u opstem slucaju). Sada postoji problem
# evaluiranje tog modela, jer vise nemamo na raspalaganju podatke koje on nije
# video. Na trening se model trenirao. Iako na test skupu nije bilo eksplicitnog
# prilagodjavanja modela podacima, ipak smo ga koristili za uporedjivanje modela
# i odabir k, pa smo implicitno odabrali najbolje k koje odgovara test podacima.
# Za evaluaciju u ovakvom slucaju se vrsi dodatna podela, izdvaja se validacioni
# skup.


## Ceo proces rada pri koriscenju validacionog skupa:
# Zelimo da odaberemo bolji model izmedju 3 modela, i da ga nakon toga evaluiramo
# Modele treniramo na trening skupu, uporedjujemo na validacionom.
# Nakon sto izaberemo najbolji, treniramo ga na trening+validacionom, evaluiramo
# na testu.
# Modeli koje koristimo su logisticka regresija, knn sa 5 suseda, i knn sa 11 suseda.
# Modele cemo da poredimo na osnovu tacnosti klasifikacije.
# Podaci iz baze Sonar (biblioteka mlbench) sadrze 60 numerickih prediktora prikupljenih
# sonarom. Ciljna promenljiva je podatak da li je skenirani objekat kamen ili mineral.
data(Sonar)
View(Sonar) #podaci su poprilicno neinterpretabilni
N = nrow(Sonar) #208
p = ncol(Sonar) - 1 #60
summary(Sonar$Class) #relativno su ujednacene klase
# Posto imamo mnogo atributa na relativno malo podataka, postoji realna opasnost
# od preprigodjavanja
library(corrplot)
corrplot(cor(Sonar[,-61]))
# Vidimo jake korelacije izmedju bliskih podataka, kao i oblasti jakih negativnih korelacija

data = Sonar
# Podelicemo podatke na 60% trening, 20% validacioni i 20% test
validation_size = 0.2
test_size = 0.2

# Postavljamo seed da bi eksperiment mogao da se reprodukuje
set.seed(217)
# Delimo podatke na trening, validacioni i test skup
# Prvo delimo na trening i ostatak
train_idx = createDataPartition(data$Class, p = (1- validation_size - test_size), list = FALSE)
train <- data[train_idx,]
rest_of_data = data[-train_idx,]

# Ostatak delimo na validacioni i test
val_idx = createDataPartition(rest_of_data$Class, p = validation_size/(validation_size + test_size), list = FALSE)
val = rest_of_data[val_idx,]
test = rest_of_data[-val_idx,] 

# Razdvajanje prediktora i ciljne promenljive
X_train = train[,-61]
y_train = as.numeric(train[, 61]) - 1

X_val = val[,-61]
y_val = as.numeric(val[, 61]) - 1

X_test = test[,-61]
y_test = as.numeric(test[, 61]) - 1

# Prediktori su napravljeni tako da budu u rasponu od 0 do 1, nema potrebe za skaliranjem


# Treniranje na treningu i uporedjivanje tri modela na validacionom skupu
log_model = glm(y_train ~., data = X_train, family = 'binomial')
# Dobijamo warning da algoritam nije iskonvergirao
# Takodje je prijavljeno da je model predvideo verovatnoce jednake 0 ili 1
# Razlog za obe pojave je cesto da logisticka regresija moze savrseno da 
# razdvoji podatke. Posto u praksi realna slika nije da su podaci savrseno 
# razdvojiji, ovo je obicno znak preprilagodjavanja.
# Regularizacione metode za linearne modele koje cemo kasnije da radimo mogu 
# ovakav problem da rese
y_pred_train = as.numeric(log_model$fitted.values>0.5)
mean(y_pred_train == y_train)

predicted_proba = predict(log_model, newdata = X_val, type = "response")
y_pred = as.numeric(predicted_proba > 0.5)
mean(y_pred == y_val) #0.6829268
# Mozete da proverite za vezbu da je tacnost klasifikacije na trening skupu 1,
# odnosno da se model izrazito preprilagodio


predicted_proba = predict(log_model, newdata = X_val, type = "response")
y_pred = as.numeric(predicted_proba > 0.5)
mean(y_pred == y_val)




knn1_y_pred = knn(X_train, X_val, y_train, k = 5)
mean(knn1_y_pred == y_val) #0.7560976

knn2_y_pred = knn(X_train, X_val, y_train, k = 11)
mean(knn2_y_pred == y_val) #0.7317073

# Najbolji je knn sa 5 suseda
# Sada treba da ga evaluiramo nepristrasno
# Treniranje na trening + validacija, evaluiranje na test skupu
X_train_val = rbind(X_train, X_val)
y_train_val = c(y_train, y_val)

knn1_full_y_pred = knn(X_train_val, X_test, y_train_val, k = 5)
mean(knn1_full_y_pred == y_test) #0.8292683
# Dobijamo iznenadjujuce odlicne performanse na test skupu
# Moguc razlog za to je sto je model sad istreniram na 80% podataka (trening + 
# validacioni), dok je pri racunanju tacnosti na validacionom skupu model bio
# istreniram na samo 60% podataka (trening skup)