rN_mesavina = function(n, mu, sigma, eps) {
  pom = sample(1:2, n, replace = T, prob = c(1 - eps, eps))
  x = rnorm(n, mu[pom], sigma[pom])
  return(x)
}

K_Gausovo = function(t) {
  exp(- t^2 / 2) / sqrt(2 * pi)
}

K_bikvadratno = function(t) {
  ifelse(abs(t) < 1, 15 / 16 * (1 - t^2)^2, 0)
}

K_trougaono = function(t) {
  ifelse(abs(t) < 1, 1 - abs(t), 0)
}

K_Epanecnikovo = function(t) {
  ifelse(abs(t) < sqrt(5), 3 / 4 / sqrt(5) * (1 - t^2 / 5), 0)
}

K_Epanecnikovo_1 = function(t) {
  ifelse(abs(t) < 1, 3 / 4 * (1 - t^2), 0)
}

K_uniformno = function(t) {
  ifelse(abs(t) < 1, 1 / 2, 0)
}

ocena_gustine_jezgrom = function(x, podaci, h, jezgro = "Gausovo") {
  fn = rep(0, length(x))
  
  if (jezgro == "Gausovo") {
    for (j in 1:length(x)) {
      fn[j] = sum(K_Gausovo((x[j] - podaci) / h))
    }
  }
  else if (jezgro == "bikvadratno") {
    for (j in 1:length(x)) {
      fn[j] = sum(K_bikvadratno((x[j] - podaci) / h))
    }
  }
  else if (jezgro == "trougaono") {
    for (j in 1:length(x)) {
      fn[j] = sum(K_trougaono((x[j] - podaci) / h))
    }
  }
  else if (jezgro == "Epanecnikovo") {
    for (j in 1:length(x)) {
      fn[j] = sum(K_Epanecnikovo((x[j] - podaci) / h))
    }
  }
  else if (jezgro == "uniformno") {
    for (j in 1:length(x)) {
      fn[j] = sum(K_uniformno((x[j] - podaci) / h))
    }
  }
  return(fn / (n * h))
}

n = 100
par(mfrow = c(1, 2))

for (a in c(1 / 2, 2)) {
  podaci = rN_mesavina(n, mu = c(a, -a), sigma = c(1, 1), eps = 1 / 2)
  x = seq(-4, 4, length.out = 1000)
  
  # Silvermanovo pravilo
  h = 0.9 * min(sd(podaci), IQR(podaci) / 1.34) * n^(-1/5)
  
  plot(x, ocena_gustine_jezgrom(x, podaci, h, "Gausovo"), 
       type = "n", xlab = "", ylab = "", xlim = c(-4, 4),
       ylim = c(0, 0.6), main = "Ocena gustine jezgrom")
  lines(x, ocena_gustine_jezgrom(x, podaci, h, "Gausovo"), type = "l", col = 1)
  lines(x, ocena_gustine_jezgrom(x, podaci, h, "bikvadratno"), type = "l", col = 2)
  lines(x, ocena_gustine_jezgrom(x, podaci, h, "trougaono"), type = "l", col = 3)
  lines(x, ocena_gustine_jezgrom(x, podaci, h, "Epanecnikovo"), type = "l", col = 4)
  lines(x, ocena_gustine_jezgrom(x, podaci, h, "uniformno"), type = "l", col = 5)
  legend(1.2, 0.5, legend=c("Gausovo jezgro", "Bikvadratno jezgro", "Trougaono jezgro", "Epanehnikovo jezgro", "Uniformno jezgro"),  
         col = 1:5, lty = rep(1, 5), cex = 0.3, box.lty = 1)
}