x<-rbeta(100, shape1 = 3, shape2 = 2)*5
n<-length(x)
y<-seq(0, 5.5, .01) #tacke u kojima se ocenjuje gustina
r<-7
M<-2^r #broj tacaka u kojima ce se odrediti ocena gustine 
h<-0.9*min(IQR(x)/1.34,sd(x))*n^(-1/5) #Silvermanova ocena koja se koristi kao polazna vrednost
#(a,b) interval u kojem se nalaze vrednosti
a<-min(x)-3*h
b<-max(x)+3*h
delta<-(b-a)/M
t<-a+delta*0:(M-1) #podeone tacke intervala
#ako X upadne u interval [t_k,t_{k+1}] dodeljuju se tezine (t_{k+1}-X)/(n*delta^2) t_k-u i 
#(X-t_k)/(n*delta^2) t_{k+1}, a onda se saberu tezine koje odgovaraju svakom t_k za svaku tacku X_i
ksi<-rep(0,length(t)) #niz u koji ce se smestati tezine za t_k
ksi1<-rep(0,length(t)) #niz u koji ce se smestati tezine za t_{k+1}
for(k in 1:(length(t)-1)){
  for(i in 1:n){
    if(x[i]>=t[k] && x[i]<=t[k+1]){
      ksi[k]<-ksi[k]+(t[k+1]-x[i])/(n*delta^2)
      ksi1[k+1]<-ksi1[k+1]+(x[i]-t[k])/(n*delta^2)
    }
  }
}
Ksi<-ksi+ksi1 #tezine za svako t_k
#Y odredjujemo Furijeovom transformacijom; zbog nacina na koji je definisano fft u R-u (help(fft))
#ovde nam je potreban inverz fft podeljen sa M
Y<-fft(Ksi, inverse = TRUE)/M 
s<-2*pi*1:M/(b-a)

psi<-exp(-h^2*s^2/2)*Y
f_ocena<-fft(psi) #ocena gustine
plot(t,f_ocena,type="l")
#ako je potrebno oderediti ocenu za jos neko h, ponoviti odredjivanje psi i f_ocena za konkretno h

#odredjivanje funkcije M1 za h iz intervala (0.25h,1.5h); ako se minimum funkcije M1 nalazi na 
#jednoj od granica onda se siri interval na tu stranu dok se ne dobije lokalni minimum

H<-seq(0.25*h,10*h,(10*h-0.25*h)/100)
M1<-c()
for(i in 1:length(H)){
  M1[i]<-(b-a)*sum(exp(-H[i]^2*s^2)-2*exp(-H[i]^2*s^2/2)*abs(Y)^2)+2/(n*H[i]*sqrt(pi)) 
}
plot(H, M1, type = "l")
h1<-H[which.min(M1)]

f<-c()
for(j in 1:length(y)){
  suma<-0
  for(i in 1:n){
    suma<-suma+dbeta(x[i]/5,y[j]/h1+1,(1-y[j])/h1+1)
  }
  f[j]<-suma
}
fh<-f/n #ocena gustine jezgrom
plot(sort(x),dbeta(sort(x), shape1 = 3, shape2 = 2)*5,type = 'l')
lines(y, fh)