#odredjivanje optimalnog h unakrsnom proverom uz pomoc Furijeove transformacije
X<-read.table("C:/Users/Marija/Desktop/OPS/SDSS1.txt")
x<-X[,2]
n<-length(x)
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,1.5*h,(1.5*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")
H[which.min(M1)]
psi<-exp(-H[which.min(M1)]^2*s^2/2)*Y
f_ocena<-fft(psi) #ocena gustine sa optimalnom glatkoscu
plot(t,f_ocena,type="l")