source("ocena_gustine_jezgrom_pom.r")

n = 100
podaci = rexp(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)

x = seq(-0.5, 5, length.out = 100)
# standardna ocena
plot(x, ocena_gustine_jezgrom(x, podaci, h_sil), type = "l",
     xlab = "", ylab = "", main = "Ocena Gausovim jezgrom")
abline(v = 0, col = "red")

# a) refleksija podataka (c = 0)
ocena_jezgrom_sa_refleksijom = function(x, podaci, h) {
  fn = rep(0, length(x))
  
  for (j in 1:length(x)) {
    if (x[j] >= 0) {
      fn[j] = ocena_gustine_jezgrom(x[j],  podaci, h) +
        ocena_gustine_jezgrom(-x[j], podaci, h)
    } else {
      fn[j] = 0
    }
  }
  fn
}

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_sil, 2.5 * h_sil, length.out = 100)
  MMV = CV(podaci, h1)
  plot(h1, MMV, type = "l", main = "MMV sa unakrsnom proverom")
  (hopt = h1[which.max(MMV)])
}

fn_ref = ocena_jezgrom_sa_refleksijom(x, podaci, hopt)
# teorijska gustina
plot(x, dexp(x), type = "l", xlab = "", ylab = "",
     main = "Ocena Gausovim jezgrom sa refleksijom", col = 4)
lines(x, fn_ref, col = 5)
lines(x, ocena_gustine_jezgrom(x, 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(x, podaci, h) {
  fn = rep(0, length(x))
  
  for (j in 1:length(x)) {
    if (x[j] >= 0) {
      for (i in 1:length(podaci)) {
        fn[j] = fn[j] +
          dgamma(podaci[i], shape = x[j] / h + 1, rate = 1 / h) / n
      }
    } else {
      fn[j] = 0
    }
  }
  fn
}

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_sil, 0.25 * h_sil, length.out = 100)
  MMV = CV(podaci, h1, ocena_gama_jezgrom)
  plot(h1, MMV, type = "l")
  (hopt = h1[which.max(MMV)])
}

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