source("ocena_gustine_jezgrom_pom.r")

n = 100
podaci = rnorm(n)

# Silvermanova ocena
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)

konvolucija = function(K) {
  pom = function(t) {
    integrate(function(u) {K(t - u) * K(u)}, -Inf, Inf)$value
  }
  return(Vectorize(pom))
}
K2 = konvolucija(K_Gausovo)

# funkcija koju zelimo da minimizujemo
M1 = function(podaci, h) {
  M_pom = rep(0, length(h))
  
  for (i in 1:n) {
    for (j in 1:n) {
      M_pom = M_pom + K2((podaci[i] - podaci[j]) / h) -
        2 * K_Gausovo((podaci[i] - podaci[j]) / h)
    }
  }
  
  M_pom = M_pom / (n^2 * h) + 2 * K_Gausovo(0) / (n * h)
  return(M_pom)
}

Mh = M1(podaci, h)
plot(h, Mh, type = 'l', main = "")

(hopt = h[which.min(Mh)])

# ako se minimum funkcije Mh nalazi na jednoj od granica 
# onda se interval siri na tu stranu dok se ne dobije lokalni minimum
if (hopt == h[length(h)]) {
  h1 = seq(h_sil, 2 * h_sil, length.out = 100)
  Mh1 = M(podaci, h1)
  plot(h1, Mh1, type = 'l', main = "")
  (hopt = h1[which.min(Mh1)])
}

# tacke u kojima ocenjujemo gustinu
x = seq(-4, 4, length.out = 100)
plot(x, ocena_gustine_jezgrom(x, podaci, hopt), type = "l",
     xlab = "", ylab = "", xlim = c(-4, 4), main = "Ocena gustine jezgrom")