library(kedd)
dens <- function(x_mix){
  K <- c()
  d <- c()
  x_min = min(x_mix)
  x_max = max(x_mix)
  h.deri = bw.nrd0(x_mix)    
  d_deri2<-dkde(x_mix,c(seq(x_min-h.deri, x_max+h.deri, 
                            by=0.1)),h=h.deri,kernel="epanechnikov",deriv.order=2)$est.fx
  d_deri2_sq<-d_deri2^2
  tot_d_deri2<-sum(d_deri2_sq)*0.1
  h_AMISE<-((3/5)/(1/25*length(x_mix)*tot_d_deri2))^(1/5)
  n = length(x_mix)
  grid<-c(seq(x_min-h_AMISE, x_max+h_AMISE, by = 0.1))
  for (i in 1:length(grid)) {
    for (j in 1:n) {
      u = (grid[i] - x_mix[j])/(h_AMISE)
      if(abs(u) <= sqrt(5)){
        K[j] = sqrt(1/5)*3/4*(1-1/5*u^2)}
      else{K[j] = 0}}
    d[i] = sum(K)/(h_AMISE*n)}
  return(d)}
