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

# a)
n = 100
podaci = rnorm(n)

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

konvolucija = function(K) {
  pom = function(z) {
    integrate(function(y) {K(z-y)*K(y)}, -Inf, Inf)$value
  }
  return(Vectorize(pom))
}

K2 = konvolucija(K_Gausovo)

# funkcija koju zelimo da minimizujemo (na H) 
J_hat = function(x, H) {
  J = rep(0, length(H))
  
  for (i in 1:n) {
    for (j in 1:n) {
      J = J + K2((x[i] - x[j])/H) - 2*K_Gausovo((x[i] - x[j])/H)
    }
  }
  
  J = J/(n^2*H) + 2*K_Gausovo(0)/(n*H)
  return(J)
}

J = J_hat(podaci, H)
plot(H, J, type = 'l', main = "MNK sa unakrsnom proverom")
(hopt = H[which.min(J)])

# ako se minimum funkcije J 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, 2*h, (2*h - h)/100)
  J1 = J_hat(podaci, H1)
  plot(H1, J1, type = 'l', main = "MNK sa unakrsnom proverom")
  (hopt = H1[which.min(J1)])
}

y = seq(-4, 4, 0.1)
plot(y, ocena_gustine_jezgrom(y, podaci, hopt), type = "l",
     xlab = "", ylab = "", xlim = c(-4, 4), main = "Ocena gustine jezgrom")

# b)
MNK_f_transf = function(x, h, M) {
  a = min(x) - 3*h
  b = max(x) + 3*h
  
  # 1. diskretizacija podataka da bi se odredio niz tezina ksi_k
  delta = (b - a)/M
  t = a + delta*0:(M-1)
  
  # tezine dodeljene granicama intervala
  ksi = rep(0, length(t))
  for (k in 1:(length(t) - 1)) {
    for (j in 1:n) {
      if (x[j] >= t[k] && x[j] < t[k+1]) {
        ksi[k] = ksi[k] + (t[k+1] - x[j])/(n*delta^2)
        ksi[k+1] = ksi[k+1] + (x[j] - t[k])/(n*delta^2)
      }
    }
  }
  
  #2. brza Furijeova transformacija za odredjivanje Y-ona
  Y = sapply((-M/2):(M/2), function(l) {
            sum(ksi*exp(1i*2*pi*0:(M-1)*l/M))
            })/M
  return(Y)
}

h_mnk = function(x, Y, s, H, a, b) {
  Y = Y[(M/2+2):(M+1)]
  s = s[(M/2+2):(M+1)]
  J = rep(0, length(H))
  for(k in 1:length(H)){
    J[k] = 2*(b-a)*sum((exp(-H[k]^2*s^2) - 2*exp(-H[k]^2*s^2/2))
                       *Mod(Y)^2) + 2/(n*H[k]*sqrt(2*pi)) - 1
  }
  return(J)
}

h_mnk = function(x, Y, s, H, a, b) {
  Y = Y[1:(M/2)]
  s = s[1:(M/2)]
  
  J = rep(0, length(H))
  for(k in 1:length(H)){
    J[k] = 2*(b-a)*sum((exp(-H[k]^2*s^2) - 2*exp(-H[k]^2*s^2/2))*Mod(Y)^2) + 
      2/(n*H[k]*sqrt(2*pi)) - 1
  }
  return(J)
}

a = min(podaci) - 3*h
b = max(podaci) + 3*h
M = 2^7
s = 2*pi*(-M/2):(M/2)/(b-a)
J = h_mnk(x, MNK_f_transf(podaci, h, M), s, H, a, b)
plot(H, J, type = 'l', main = "MNK sa unakrsnom proverom - Furijeova transformacija")

(hopt = H[which.min(J)])

if (hopt == H[length(H)]) {
  H1 = seq(h, 2*h, (2*h - h)/100)
  J1 = h_mnk(x, MNK_f_transf(podaci, h, M), s, H1, a, b)
  plot(H1, J1, type = 'l', main = "MNK sa unakrsnom proverom - Furijeova transformacija")
  (hopt = H1[which.min(J1)])
}