source("ocena_gustine_jezgrom_pom.r")

n = 100
podaci = rnorm(n)

h_sil = 0.9 * min(sd(podaci), IQR(podaci) / 1.34) * n^(-1/5)
h = seq(0.25 * h_sil, 1.5 * h_sil, length.out = 100)

# funkcija koju zelimo da maksimizujemo
CV = function(uzorak, h, ocena_gustine = ocena_gustine_jezgrom){
  # po redovima: red i -> ocene gustine bez i-tog elementa, i = 1:n
  f_i = matrix(0, n, length(h))
  
  for (i in 1:n) {
    # x - tacke u kojima racunamo ocenu gustine, podaci - uzorak iz koga dobijamo ocenu
    f_i[i, ] = sapply(h, ocena_gustine, x = uzorak[i], podaci = uzorak[-i])
  }
  
  cv = rep(0, length(h))
  for (j in 1:length(h)) {  
    cv[j] = sum(log(f_i[, j])) / n
  }
  return(cv)
}

MMV = CV(podaci, h, ocena_gustine_jezgrom)
plot(h, MMV, type = "l", main = "")

if (hopt == h[length(h)]) {
  h1 = seq(1.5 * h_sil, 2.5 * h_sil, length.out = 100)
  MMV = CV(podaci, h1, ocena_gustine_jezgrom)
  plot(h1, MMV, type = "l", main = "")
  (hopt = h1[which.max(MMV)])
}

x = seq(-4, 4, length.out = 100)
plot(x, ocena_gustine_jezgrom(x, podaci, hopt), type = "l",
     xlab = "", ylab = "", xlim = c(-4, 4), main = "Ocena gustine jezgrom")