###### ilustracija CV metode u kNN regresiji #######

## model je Y=f_true(X) + epsilon sa E[epsilon]=0 i Var[epsilon]=sigma^2

## Primjer 1 -- f_true nije linearna + velika varijanca
f_true=function(x){return(x^2 + 5*cos(x^2))}
sigma=10

# # Primjer 2 -- f_true je linearna
# f_true=function(x){return(x)}
# sigma=2
# 
# ## Primjer 3 -- f_true nije linearna + mala varijanca
# f_true=function(x){return(x^2 + 5*cos(x^2))}
# sigma=1

library("FNN")

set.seed(5) ## kako bi mogli ponoviti rezultate simulacije
n=100
x_train=runif(n,1,4)
eps_train=rnorm(n,0,sigma)
y_train=f_true(x_train)+eps_train

plot(x_train, y_train, type = "p", xlab="x", ylab="y")
t=seq(1, 4, by=0.01)
points(t, f_true(t), type = "l")

moj_CV=function(k, K,x,y){ ## k je parametar za CV, a K za metodu najblizih susjeda
  n=length(x)
  r=n/k ## velicina blokova, pretpostavljamo da je to prirodan broj
  CV=numeric(k)
  for (j in 1:k) {
    leave_out=((j-1)*r+1):(j*r)   ## buduci da znamo da su podaci njd, mozemo blokove odabrati deterministicki
    tr=x[-leave_out] 
    te=x[leave_out]
    
    y_tr=y[-leave_out]
    y_te=y[leave_out]
    ## (tr, y_tr) predstavljaju \tau^{-j}
    ## (te,y_te) predstavljaju \tau_j koji igra ulogu skupa za validaciju
    knn_pred=knn.reg(train=tr, y=y_tr, test = cbind(te), k=K)$pred
    CV[j]=mean((y_te-knn_pred)^2)
  }
  SE = sd(CV)/sqrt(k) ## (naivni) procjenitelj standardne greske 
  return(c(mean(CV), SE))
}

k=10
r=n/k ## velicina bloka
## moci cemo napraviti CV procjenu za KNN metodu samo u slucaju kada je
## K<= n-r jer je (n-r) duljina trening skupa u svakom koraku CV algoritma
CV=numeric(n-r-1)
SE=numeric(n-r-1)

## racunamo Cv procjenu greske i procjenu njene standardne greske za svaki K
for (K in 1:(n-r-1)) {
  z=moj_CV(k,K,x_train, y_train)
  CV[K]=z[1]
  SE[K]=z[2]
}

## testna greska L(f^hat_K)
N=10000
x_test=runif(N,1,4)
eps_test=rnorm(N,0,sigma)
y_test=f_true(x_test)+eps_test 
test_error_knn=numeric(n-r-1)
for (K in 1:(n-r-1)) {
  knn_fit=knn.reg(train=x_train, test=cbind(x_test),  y=y_train, k=K)
  test_error_knn[K]=mean((y_test-knn_fit$pred)^2)
}

### crtanje
M=max(test_error_knn, CV)
m=min(test_error_knn, CV)

Kna_min_1=log(((n-r-1):1)^(-1))
plot(Kna_min_1, rev(test_error_knn), type="b",pch=16, cex=0.5, ylim=c(m,M),
     col="orange", xlab="log(1/K) (fleksibilnost kNN modela)", ylab="")
knn_min=which.min(test_error_knn)
points(-log(knn_min), min(test_error_knn), pch=15, col="orange")

points(Kna_min_1, rev(CV), type="b",pch=16, cex=0.5, col="blue")
CV_min=which.min(CV) ## 14
points(-log(CV_min), min(CV), pch=15, col="blue")
abline(h=min(CV)+SE[CV_min],lty=2)  ## min(CV) + jedna standardna greska

CV_1SE_rule_ind=max(which(CV<=min(CV) + SE[CV_min]))
## trazimo najveci K (dakle najmanje kompleksan model) cija CV procjena
## greske ne prelazi minimalnu + SE 
CV_1SE_rule_ind  ## 29
CV_1SE_rule=CV[CV_1SE_rule_ind]

points(-log(CV_1SE_rule_ind), CV_1SE_rule, pch=15, col="lightblue")

legend(Kna_min_1[30], M, legend=c("Test error", "CV error"), 
       col=c("orange", "blue"),
       lty=c(1,1))

# ## ireducibilna greska
# abline(h=sigma^2, col="red", lty=4)

## slucaj 1, k=10 
## -> oblik CV krivulje slican obliku krivulje stvarne greske 
## (jedino bitno ako nas zanima samo odabir modela)


## isprobati druge primjere i parametre n,k