##############################################################
######## MIXTURE OF NORMAL DENSITIES #########################
##############################################################

set.seed(0)

library(nor1mix)

#-- Samples

mas_1<-rnorMix(100,MW.nm1)
mas_2<-rnorMix(100,MW.nm2)
mas_3<-rnorMix(100,MW.nm3)
mas_4<-rnorMix(100,MW.nm4)
mas_5<-rnorMix(100,MW.nm5)
mas_6<-rnorMix(100,MW.nm6)
mas_7<-rnorMix(100,MW.nm7)
mas_8<-rnorMix(100,MW.nm8)
mas_9<-rnorMix(100,MW.nm9)
mas_10<-rnorMix(100,MW.nm10)
mas_11<-rnorMix(100,MW.nm11)
mas_12<-rnorMix(100,MW.nm12)
mas_13<-rnorMix(100,MW.nm13)
mas_14<-rnorMix(100,MW.nm14)
mas_15<-rnorMix(100,MW.nm15)

mas<-t(cbind(sample_1,sample_2,sample_3,sample_4,sample_5,sample_6,sample_7,sample_8,sample_9,sample_10,sample_11,sample_12,sample_13,sample_14,sample_15))


#-- Theorical densities

density_1<-dnorMix(mas_1,MW.nm1)
density_2<-dnorMix(mas_2,MW.nm2)
density_3<-dnorMix(mas_3,MW.nm3)
density_4<-dnorMix(mas_4,MW.nm4)
density_5<-dnorMix(mas_5,MW.nm5)
density_6<-dnorMix(mas_6,MW.nm6)
density_7<-dnorMix(mas_7,MW.nm7)
density_8<-dnorMix(mas_8,MW.nm8)
density_9<-dnorMix(mas_9,MW.nm9)
density_10<-dnorMix(mas_10,MW.nm10)
density_11<-dnorMix(mas_11,MW.nm11)
density_12<-dnorMix(mas_12,MW.nm12)
density_13<-dnorMix(mas_13,MW.nm13)
density_14<-dnorMix(mas_14,MW.nm14)
density_15<-dnorMix(mas_15,MW.nm15)

densities<-t(cbind(density_1,density_2,density_3,density_4,density_5,density_6,density_7,density_8,density_9,density_10,density_11,density_12,density_13,density_14,density_15))


##############################################################
######## AUXILIAR FUNCTIONS ##################################
##############################################################

#-- KERNEL

kernel<-function(x){
return(1/sqrt(2*pi)*exp(-x^2/2))
}

#-- KERNEL ESTIMATOR

kernel_estimator<-function(x,bandwidth,mas){
return(1/(length(mas)*bandwidth)*apply(kernel(outer(x,mas,"-")/bandwidth),1,sum))
}

#-- ISE INTEGRAND

ISE_integrand<-function(x,x_points){
return((x-dnorMix(x_points,MW.nm1))^2)
}

#-- LCV

LCV<-function(bandwidth,mas){
return(sum(dnorm(outer(mas,mas,"-")/bandwidth,sd=sqrt(2)))/(length(mas)^2*bandwidth)- (4/(length(mas)*(length(mas)-1)*bandwidth)) * sum(kernel((outer(mas,mas,"-")/bandwidth)[upper.tri(matrix(0,nr=length(mas),nc=length(mas)))])))
}



#################################################################################################
######## PARAMETERS #############################################################################
#################################################################################################


B<-500

ISE_cv<-numeric(length=B)
ISE_normal<-numeric(length=B)
cross_validation_bandwidth<-numeric(length=B)
normal_scale_bandwidth<-numeric(length=B)


for(b in 1:B){


##############################################################
######### CROSS-VALIDATION BANDWIDTH #########################
##############################################################


mas<-rnorMix(100,MW.nm1)

#LCV<-function(bandwidth,mas){

#evaluating_points<-as.numeric(sort(mas)%*%(cbind(rbind(diag(1/2,length(mas)-1),rep(0,length(mas)-1)),rep(0,length(mas)))+ cbind(rbind(rep(0,length(mas)-1),diag(1/2,length(mas)-1)),rep(0,length(mas)))))[-length(mas)]
#intervals_mas<-as.numeric(sort(mas)%*%(cbind(rbind(diag(-1,length(mas)-1),rep(0,length(mas)-1)),rep(0,length(mas)))+ cbind(rbind(rep(0,length(mas)-1),diag(1,length(mas)-1)),rep(0,length(mas)))))[-length(mas)]
#return( (1/(length(mas)^2 * bandwidth^2)) * sum((apply(kernel(outer(evaluating_points,mas,"-")/bandwidth),1,sum))^2*intervals_mas)
#- (4/(length(mas)*(length(mas)-1)*bandwidth)) * sum(kernel((outer(mas,mas,"-")/bandwidth)[upper.tri(matrix(0,nr=100,nc=100))])))
#}


#CV<-numeric()
#h<-seq(0,max(dist(mas)),by=0.01)
#for (i in 1:length(h)){CV[i]= LCV(bandwidth=h[i],mas)}
#plot(h,CV,type="l")

cross_validation_bandwidth[b]<-optimize(LCV,interval=c(0,max(dist(mas))),mas=mas)$minimum


##############################################################
######### NORMAL BANDWIDTH ###################################
##############################################################


normal_scale_bandwidth[b]<-1.06*min(sd(mas),IQR(mas)/(qnorm(0.75)-qnorm(0.25)))*length(mas)^(-1/5)


##############################################################
######### ISE ################################################
##############################################################


x<-seq(min(mas),max(mas),by=0.01)

#kernel_estimator_cv<-1/(length(mas)*cross_validation_bandwidth[b])*apply(kernel((outer(x,mas,"-"))/cross_validation_bandwidth[b]),1,sum)
#kernel_estimator_normal<-1/(length(mas)*normal_scale_bandwidth[b])*apply(kernel((outer(x,mas,"-"))/normal_scale_bandwidth[b]),1,sum)

#plot(kernel_estimator_cv~x,type="l")
#curve(dnorMix(x,MW.nm1),lty=3,add=T)
#windows()
#plot(kernel_estimator_normal~x,type="l")
#curve(dnorMix(x,MW.nm1),lty=3,add=T)

ISE_cv[b] <- cross_validation_bandwidth[b]/2 * (ISE_integrand(kernel_estimator(x[1],cross_validation_bandwidth[b],mas),x[1])+ ISE_integrand(kernel_estimator(x[length(x)],cross_validation_bandwidth[b],mas),x[length(x)]))+
cross_validation_bandwidth[b] * sum(ISE_integrand(kernel_estimator(x[2:(length(x)-1)],cross_validation_bandwidth[b],mas),x[2:(length(x)-1)]))

ISE_normal[b] <- normal_scale_bandwidth[b]/2 * (ISE_integrand(kernel_estimator(x[1],normal_scale_bandwidth[b],mas),x[1])+ ISE_integrand(kernel_estimator(x[length(x)],normal_scale_bandwidth[b],mas),x[length(x)]))+
normal_scale_bandwidth[b] * sum(ISE_integrand(kernel_estimator(x[2:(length(x)-1)],normal_scale_bandwidth[b],mas),x[2:(length(x)-1)]))


}

