###### 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]
    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)
  return(c(mean(CV), SE))
}

k=10
r=n/k
CV=numeric(n-r-1) ## za K >=n-r ne mozemo provesti CV
SE=numeric(n-r-1)

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
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="l", 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="l", col="blue")
CV_min=which.min(CV)
points(-log(CV_min), min(CV), pch=15, col="blue")
abline(h=min(CV)+SE[CV_min],lty=2)  ## min(CV) + jedna standardna greska

## 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),
## ali daje i "solidnu" procjenu testne greske
## One SE rule ovdjte takoder daje skoro pa optimalan model sto se tice predikcije

## isprobati druge primjere i parametre n,k