## model je Y=f_true(X) + epsilon sa epsilon~N(0,sigma^2), i 
## epsilon nezavisna od X

## 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


## jedna relizacija skupa za ucenje duljine n=100
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")

# pravac (tj. f_hat) dobiven metodom najmanjih kvadrata
lin_fit=lm(y_train~x_train)
abline(lin_fit, col="blue")

library(FNN)
?knn.reg
K1=18 ## u Primjeru 1 ovo je k za koji se postize najmanja testna greska 
knn_fit=knn.reg(train=x_train, test = cbind(t),  y=y_train, k=K1)
points(t, knn_fit$pred, type="l", col="red")
## knn_fit$pred daje vektor predikcija za skup test

K2=2 ## 
knn_fit=knn.reg(train=x_train, test = cbind(t),  y=y_train, k=K2)
points(t, knn_fit$pred, type="l", col="darkred") # -> overfitting u Primjerima 1 i 2!

### probati i druge k-ove!

## dodavanje opisa na graf
legend(1, max(y_train), legend=c("Prava f", "Lin. model",
                                 paste("kNN (k=", K1, sep="", ")" ),  paste("kNN (k=", K2, sep="", ")" )), 
       col=c(1, "blue","red", "darkred", "red"), lty=c(1,1,1,1,1))



##### greska na skupu za ucenje ("training error") 
train_error_lin=mean(lin_fit$residuals^2)

train_error_knn=numeric(n-1)
for (K in 1:(n-1)) {
  knn_fit=knn.reg(train=x_train, test=cbind(x_train),  y=y_train, k=K)
  train_error_knn[K]=mean((y_train-knn_fit$pred)^2)
}


##### greska na testnom skupu
m=10000
x_test=runif(m,1,4)
eps_test=rnorm(m,0,sigma)
y_test=f_true(x_test)+eps_test 
## u praksi ne mozemo simulirati testni skup jer ne znamo
## distribuciju od (X,Y)

y_lin_test=predict(lin_fit, newdata = data.frame(x_train=x_test))
test_error_lin=mean((y_test-y_lin_test)^2)

test_error_knn=numeric(n-1)
for (K in 1:(n-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)
}


M=max(test_error_knn)
m=0

Kna_min_1=log(((n-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")

## greska na testnom skupu za linearni model
abline(h=test_error_lin, col="orange", lty=2)
## ireducibilna greska
abline(h=sigma^2, col="red", lty=4)


### dodavanje i greske na skupu za ucenje
points(Kna_min_1, rev(train_error_knn), type="l", col="blue")
abline(h=train_error_lin, col="blue", lty=2)

## dodavanje opisa na graf
legend(Kna_min_1[30], M, legend=c("Test error (kNN)", "Test error (Lin. model)",
                                  "Training error (kNN)", "Training error (Lin. model)", 
                                  "Irreducible error"), 
       col=c("orange", "orange","blue", "blue", "red"),
       lty=c(1,2,1,2,4))