source("96.r")
# zbog funkcije ocena_gustine_jezgrom

n = 100
podaci = rnorm(n)

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

# funkcija koju zelimo da maksimizujemo (na H)
CV = function(uzorak, H, ocena = ocena_gustine_jezgrom){
  f_i = matrix(0, n, length(H))
  
  # ocene gustine bez i-tog elementa (smestene po redovima)
  for (i in 1:n) {
    # y - tacke u kojima racunamo ocenu gustine, x - uzorak iz koga dobijamo ocenu
    f_i[i, ] = sapply(H, ocena, y = uzorak[i], x = 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 = "MMV sa unakrsnom proverom")
(hopt = H[which.max(MMV)])

if (hopt == H[length(H)]) {
  H1 = seq(1.5*h, 2.5*h, (2.5*h-1.5*h)/100)
  MMV = CV(podaci, H1, ocena_gustine_jezgrom)
  plot(H1, MMV, type = "l", main = "MMV sa unakrsnom proverom")
  (hopt = H1[which.max(MMV)])
}

y = seq(-4, 4, 0.01)
plot(y, ocena_gustine_jezgrom(y, podaci, hopt), type = "l", 
     xlab = "", ylab = "", main = "Ocena gustine jezgrom")