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

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

y = seq(-1, 5, 0.01)
# standardna ocena
plot(y, ocena_gustine_jezgrom(y, podaci, h), type = "l",
     xlab = "", ylab = "", main="Ocena normalnim jezgrom")
abline(v = 0, col = 'red')

# a) refleksija podataka (c = 0)
ocena_jezgrom_sa_refleksijom = function(y, x, h) {
  f = rep(0,length(y))
  
  for(j in 1:length(y)) {
    f[j] = sum(sapply((y[j] - x)/h, K_Gausovo) + 
                 sapply((-y[j] - x)/h, K_Gausovo))
  }
  f = f/(n*h)
  return(f)
}

MMV = CV(podaci, H, ocena_jezgrom_sa_refleksijom)
plot(H, MMV, type = "l", main = "MMV sa unakrsnom proverom")
(hopt = H[which.max(MMV)])

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

# da bi bila gustina, treba za negativne y da bude ocena 0, a inace dobijena vrednost
f_ref = c(rep(0, length(which(y < 0))), 
        ocena_jezgrom_sa_refleksijom(y[-which(y < 0)], podaci, hopt))
# teorijska gustina
plot(y, dexp(sort(y)), type = "l", xlab = "", ylab = "",
     main = "Ocena Gausovim jezgrom sa refleksijom", col = 4)
lines(y, f_ref, col = 5)
lines(y, ocena_gustine_jezgrom(y, podaci, hopt), col = 6)
abline(v = 0, col = 2)
legend(3, 0.5, legend = c("Teorijska raspodela", "Ocena sa refleksijom", "Ocena Gausovim jezgrom"),  
       col = 4:6,lty = rep(1,5), cex = 0.5, box.lty = 1)

# b)
ocena_gama_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] + dgamma(x[i], y[j]/H + 1, rate = 1/H)/n 
    }
  }
  return(f)
}

MMV = CV(podaci, H, ocena_gama_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_gama_jezgrom)
  plot(H1, MMV, type = "l")
  (hopt = H1[which.max(MMV)])
}

f_gamma = c(rep(0, length(which(y < 0))),
                        ocena_gama_jezgrom(y[-which(y < 0)], podaci, hopt))
plot(y, dexp(y), type = 'l', xlim = c(-1, 5), ylim = c(0, 1.1),
     xlab = "", ylab = "", main = "Ocena gama jezgrom", col = 4)
lines(y, f_gamma, col = 5)
lines(y, ocena_gustine_jezgrom(y, podaci, H[which.max(CV(podaci, H))]), col = 6)
abline(v = 0, col = 2)
legend(3, 0.5, legend = c("Teorijska raspodela", "Ocena sa refleksijom", "Ocena Gausovim jezgrom"),  
       col = 4:6,lty = rep(1,5), cex = 0.5, box.lty = 1)