#odredjivanje optimalnog h metodom najmanjeg kvadrata sa unakrsnom proverom ako se koristi Gausovo jezgro
x<-rnorm(100,0,1)
n<-length(x)
Kn<-function(t){
  return(exp(-t^2/2)/sqrt(2*pi))
}
h<-0.9*min(IQR(x)/1.34,sd(x))*n^(-1/5)
#konvolucija dve funkcije
dSumXY <- function(dX, dY){
  #konvolucija raspodela
  dReturn <- function(z){
    integrate( function(x,z){
      dX(x) * dY(z - x) },
      -Inf, Inf, z)$value
  }
  #vraca vektor
  dReturn <- Vectorize(dReturn)
  return(dReturn)  
}
convXY<-dSumXY(Kn, Kn)
#preporuceni interval na kojem se odredjuje minimum funkcije M1: (0.25h,1.5h); 
#odredjuje se vrednost funkcije za proizvoljan broj tacaka (ovde je uradjeno za 101)
H<-seq(0.25*h,1.5*h,(1.5*h-0.25*h)/100)
M<-rep(0,length(H))
for(i in 1:n){
  for(j in 1:n){
    z<-(x[i]-x[j])/H
    M<-M+convXY(z)-2*Kn(z)
  }
}
#M1<-n^(-2)*h^(-1)*sum(K^*((Xi-Xj)/h))+2*K(0), gde je K^*=K^(2)-2K, gde je K^(2) konvolucija jezgra sa samim sobom
M1<-M/(n^2*H)+2*Kn(0)/(n*H)
plot(H, M1, type = "l")
H[which.min(M1)]
#ako se minimum funkcije M1 nalazi na jednoj od granica onda se siri interval na tu stranu dok se 
#ne dobije lokalni minimum