source("ocena_gustine_jezgrom_pom.r")

n = 100
podaci = rbeta(n, 3, 2) * 5

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

ocena_beta_jezgrom = function(x, podaci, h) {
  fn = rep(0, length(x))
  
  for (j in 1:length(x)) {
    if(x[j] < 0 || x[j] > 5) {
      fn[j] = 0
    }
    else {
      for (i in 1:length(podaci)) {
        fn[j] = fn[j] + dbeta(podaci[i] / 5, x[j] / (5 * h), (5 - x[j]) / (5 * h)) / (n * 5) 
      }
    }
  }
  return(fn)
}

MMV = CV(podaci, h, ocena_beta_jezgrom)
plot(h, MMV, type = "l", main = "MMV sa unakrsnom proverom")

(hopt = h[which.max(MMV)])

if(hopt == h[1]) {
  h1 = seq(0.05 * h_sil, 0.25 * h_sil, length.out = 100)
  MMV = CV(podaci, h1, ocena_beta_jezgrom)
  plot(h1, MMV, type = "l", main = "MMV sa unakrsnom proverom")
  (hopt = h1[which.max(MMV)])
}

fn_beta = ocena_beta_jezgrom(x, podaci, hopt)
# teorijska gustina
plot(x, sapply(x / 5, dbeta, 3, 2) / 5, type = 'l',
     xlab = "", ylab = "", main = "Ocena beta jezgrom", col = 1, ylim = c(0, 0.5))
lines(x, fn_beta, col = 2)
lines(x, ocena_gustine_jezgrom(x, podaci, h[which.max(CV(podaci, h))], jezgro = "Gausovo"), col = 3)
abline(v = 0, col = 4)
legend(-0.3, 0.3, legend = c("Teorijska raspodela", "Ocena beta jezgrom", "Ocena Gausovim jezgrom"),  
       col = 1:3, lty = rep(1, 5), cex = 0.5, box.lty = 1)