#### Primjer 2.4 s predavanja

library("FNN")

set.seed(12)
n=100
p=5000
x_train=matrix(rnorm(n*p), nrow=n)
y_train=rbinom(n,1,1/2)

korelacije=numeric(p)

for(i in 1:p){
  korelacije[i]=cor(x_train[,i],y_train)
}

plot(korelacije, type="p", pch=16,cex=0.5)
## za svaki i, stvarna korelacija je 0
## ipak uzoracka korelacije je slucajna varijabla te ce za neke od kovarijata
## zbog slucajnosti ta korelacija biti relativno velika (pogotovo ako je p velik
## kao sto je ovdje slucaj)

## trazimo m=100 kovarijata s najvecom uzorackom korelacijom sa y
ind=sort(korelacije, index.return=TRUE, decreasing=TRUE)$ix
m=100
plot(korelacije[ind[1:m]], type="p", pch=16,cex=0.5)

## zadrzavamo samo tih m kovarijata
x_novi=x_train[,ind[1:m]]

## zanima nas greska 1NN metode za klasifikaciju
## kasnije cemo to raditi detaljno, ali ova metoda je vrlo jednostavna:
## za svaki x iz R^m kao predikciju uzmi klasu (0 ili 1) najblizeg susjeda
## od x u x_novi.


## 5-fold CV procjena greske -- krivi nacin
r=n/5
CV=numeric(5)
for (j in 1:5) {
  leave_out=((j-1)*r+1):(j*r)   
  ## buduci da znamo da su podaci njd, mozemo blokove odabrati deterministicki
  tr=x_novi[-leave_out,]
  te=x_novi[leave_out,]
  y_tr=y_train[-leave_out]
  y_te=y_train[leave_out]
  knn_pred=knn(train=tr, cl=y_tr, test = te, k=5)
  CV[j]=mean(knn_pred!=y_te) ## 0-1 gubitak -- jednostavno brojimo prosjecan
                              ## broj krivih predikcija
}
CV_error=mean(CV)
CV_error  ## 0.07 ?!

## problem je taj sto smo kovarijate u x_novi izabrali na temelju gledanja
## cijelog skupa za trening, odnosno svih odziva y_train!
## zbog toga, u svakom koraku CV algoritma, izostavljeni skup (te,y_te) 
## vise ne predstavlja dobar skup za validaciju (tj. nezavisan testni skup)!

## proces odabira m kovarijata na temelju uzoracke korelacije s odzivom je
## takoder dio algoritma ucenja te zato treba biti napravljen ponovno
## u svakom koraku CV algoritma!
  
## 5-fold CV -- ispravan nacin
r=n/5
CV=numeric(5)
for (j in 1:5) {
  leave_out=((j-1)*r+1):(j*r)   
  tr=x_train[-leave_out,]
  te=x_train[leave_out,]
  y_tr=y_train[-leave_out]
  y_te=y_train[leave_out]
    
  ## racunamo korelacije, ali samo na (tr,y_tr)!
  korelacije_cv=numeric(p)
  for(i in 1:p){
    korelacije_cv[i]=cor(tr[,i],y_tr)
  }
    
  ind_cv=sort(korelacije_cv, index.return=TRUE, decreasing=TRUE)$ix
    
  tr=tr[,ind_cv[1:m]]
  te=te[,ind_cv[1:m]]
    
  knn_pred=knn(train=tr, cl=y_tr, test = te, k=5)
  CV[j]=mean(knn_pred!=y_te)
}
CV_error=mean(CV)
CV_error ## 0.49


## simulacijska studija
# ponovit cemo gornji postupak M puta te na taj nacin dobiti uzoraka duljine
## M za CV provjenitelj greske u oba slucaja (algoritam traje malo duze) 
## -- pogledaj gotovu sliku CV_krivo_i_ispravno.png
M=500
CV_krivo=numeric(M)
CV_ispravno=numeric(M)

for (t in 1:M) {
  x_train=matrix(rnorm(n*p), nrow=n)
  y_train=rbinom(n,1,1/2)
  
  korelacije=numeric(p)
  
  for(i in 1:p){
    korelacije[i]=cor(x_train[,i],y_train)
  }
  
  ind=sort(korelacije, index.return=TRUE, decreasing=TRUE)$ix
  x_novi=x_train[,ind[1:m]]
  
  ## 5-fold CV procjena greske -- krivi nacin
  r=n/5
  CV=numeric(5)
  for (j in 1:5) {
    leave_out=((j-1)*r+1):(j*r)   
    ## buduci da znamo da su podaci njd, mozemo blokove odabrati deterministicki
    tr=x_novi[-leave_out,]
    te=x_novi[leave_out,]
    y_tr=y_train[-leave_out]
    y_te=y_train[leave_out]
    knn_pred=knn(train=tr, cl=y_tr, test = te, k=5)
    CV[j]=mean(knn_pred!=y_te)
  }
  CV_krivo[t]=mean(CV)
  
  ## 5-fold CV -- ispravan nacin
  r=n/5
  CV=numeric(5)
  for (j in 1:5) {
    leave_out=((j-1)*r+1):(j*r)   
    ## buduci da znamo da su podaci njd, mozemo blokove odabrati deterministicki
    tr=x_train[-leave_out,]
    te=x_train[leave_out,]
    y_tr=y_train[-leave_out]
    y_te=y_train[leave_out]
    
    ## racunamo korelacije, ali samo na (tr,y_tr)!
    korelacije_cv=numeric(p)
    for(i in 1:p){
      korelacije_cv[i]=cor(tr[,i],y_tr)
    }
    
    ind_cv=sort(korelacije_cv, index.return=TRUE, decreasing=TRUE)$ix
    
    tr=tr[,ind_cv[1:m]]
    te=te[,ind_cv[1:m]]
    
    knn_pred=knn(train=tr, cl=y_tr, test = te, k=5)
    CV[j]=mean(knn_pred!=y_te)
  }
  CV_ispravno[t]=mean(CV)
}

hist(CV_krivo, xlim=c(0,0.8), col="darkred",
     main="Primjer 2.4", xlab="CV procjena greske (crveno je krivi nacin)")
hist(CV_ispravno, add=TRUE, col="darkgreen", breaks=15)
abline(v=0.5, lwd=2)
mean(CV_ispravno) ## [1] 0.50314

## U ovom slucaju se cini da je (ispravna) CV procjena greske nepristran procjenitelj
## za ocekivanu testnu gresku (=0.5), ali i za pravu testnu gresku jer ona za svaki uzorak
## takoder jednaka 0.5