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

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

y = seq(-0.5, 5.5, 0.05) 
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)

K_beta = function(x, a, b) {
  if(a < 0 || b < 0) {
    return(0)
  }
  else {
    return(dbeta(x, a, b))
  }
}

ocena_beta_jezgrom = function(y, x, H) {
  f = rep(0,length(y))
  
  for (j in 1:length(y)) {
    for (i in 1:length(x)) {
      f[j] = f[j] + K_beta(x[i]/5, y[j]/(5*H), (5-y[j])/(5*H))/(n*5) 
    }
  }
  return(f)
}

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, 0.25*h, (0.25*h-0.05*h)/100)
  MMV = CV(podaci, H1, ocena_beta_jezgrom)
  plot(H1, MMV, type = "l", main = "MMV sa unakrsnom proverom")
  (hopt = H1[which.max(MMV)])
}

f_beta = ocena_beta_jezgrom(y, podaci, hopt)
# teorijska gustina
plot(sort(y), sapply(y/5, dbeta, 3, 2)/5, type = 'l',
     xlab = "", ylab = "", main = "Ocena beta jezgrom", col = 4)
lines(y, f_beta, col = 5)
lines(y, ocena_gustine_jezgrom(y, podaci, hopt), col = 6)
legend(-0.3, 0.35, legend = c("Teorijska raspodela", "Ocena beta jezgrom", "Ocena Gausovim jezgrom"),  
       col = 4:6,lty = rep(1,5), cex = 0.5, box.lty = 1)