##############################################################
######## 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(mas_1,mas_2,mas_3,mas_4,mas_5,mas_6,mas_7,mas_8,mas_9,mas_10,mas_11,mas_12,mas_13,mas_14,mas_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))


##############################################################
######## NORMAL KERNEL #######################################
##############################################################


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

normalkernel<-function(x,mas,bandwidth){
    nx=length(x)
    fx=numeric(nx)
    
    for (i in 1:nx){
    aux<-mean(kernel((mas-x[i])/bandwidth))
    fx[i]=aux/bandwidth
    }
    
    return(fx)
}




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


#mas<-mas_1

#integral<-function(x,bandwidth,mas){
#evaluation_kernel<-kernel(outer(x,mas,"-")/bandwidth)
#return(apply(apply(evaluation_kernel,1,function(x) outer(x,x)),2,sum))
#}

#integral<-function(x,bandwidth,mas){
#evaluation_kernel<-kernel(outer(x,mas,"-")/bandwidth)
#return(apply(convolve(evaluation_kernel,evaluation_kernel),1,sum))
#}


#return(apply(apply(apply(apply(,2,function(x){x<-matrix(x,byrow=T,nrow=mas_size)}),1,sum),nrow=mas_size),2,sum))
#evaluation_kernel[upper.tri(evaluation_kernel,diag=T)]<-0
#return(2*sum(apply(outer(kernel(outer(x,mas,"-")/bandwidth),evaluation_kernel),1,sum)))
#evaluation_kernel<-kernel(outer(x,mas,"-")/bandwidth)
#sum(apply(apply(evaluation_kernel,1,function(x) outer(x,x)),1,sum))


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)]
evaluation_kernel<-kernel(outer(evaluating_points,mas,"-")/bandwidth)
return(1/(length(mas)^2 * bandwidth^2) * sum(apply(apply(evaluation_kernel,1,function(x) outer(x,x)),2,sum)*intervals_mas)
-4/(length(mas)*(length(mas)-1)*bandwidth)*sum(kernel((outer(mas,mas,"-")/bandwidth)[upper.tri(outer(mas,mas,"-")/bandwidth,diag=F)])))

}

#cross_validation_bandwidth<-optimize(LCV,interval=c(0,max(dist(mas))),mas=mas)


##############################################################
######### NORMAL SCALE BANDWIDTH #############################
##############################################################

#normal_scale_bandwidth<-1.06*min(sqrt(var(mas)),IQR(mas)/(qnorm(0.75)-qnorm(0.25)))*length(mas)^(-1/5)


##############################################################
######## INTEGRATED SQUARED ERROR#############################
##############################################################
#kernel_estimator<-1/(length(mas)*bandwidth)*apply(kernel((outer(x,mas,-))/bandwidth),1,sum)
#para mas_1
x<-seq(min(mas_1),max(mas_1),by=0.01)

#par(mfrow=c(1,2))
#plot(normalkernel(x,mas_1,cross_validation_bandwidth$minimum))
#plot(normalkernel(x,mas_1,normal_scale_bandwidth))
M<-50
B<-50
aux1<-numeric(M)
aux2<-numeric(M)
ISE<-matrix(0,2,B)
for(b in 1:B){
    x0<-sample(x,length(mas_1),replace=T)
    cross_validation_bandwidth<-optimize(LCV,interval=c(0,max(dist(mas_1))),mas=mas_1)
    normal_scale_bandwidth<-1.06*min(sqrt(var(mas_1)),IQR(mas_1)/(qnorm(0.75)-qnorm(0.25)))*length(mas_1)^(-1/5)
    density_CV<-normalkernel(x0,mas_1,cross_validation_bandwidth$minimum)
    density_normal<-normalkernel(x0,mas_1,normal_scale_bandwidth)
    density_theorical=densities[1,]
    for(m in 1:M){
    aux1[m]<-sum((density_CV-density_theorical)^2)
    aux2[m]<-sum((density_normal-density_theorical)^2)
    }
    ISE[1,b]<- mean(aux1)
    ISE[2,b]<-mean(aux2)
    }
    

##############################################################

B<-500
M<-50
aux1<-numeric(M)
aux2<-numeric(M)
ISE_1<-matrix(0,2,B)
ISE_2<-matrix(0,2,B)
ISE_3<-matrix(0,2,B)
ISE_4<-matrix(0,2,B)
ISE_5<-matrix(0,2,B)
ISE_6<-matrix(0,2,B)
ISE_7<-matrix(0,2,B)
ISE_8<-matrix(0,2,B)
ISE_9<-matrix(0,2,B)
ISE_10<-matrix(0,2,B)
ISE_11<-matrix(0,2,B)
ISE_12<-matrix(0,2,B)
ISE_13<-matrix(0,2,B)
ISE_14<-matrix(0,2,B)
ISE_15<-matrix(0,2,B)
ISE<-t(rbind(ISE_1,ISE_2,ISE_3,ISE_4,ISE_5,ISE_6,ISE_7,ISE_8,ISE_9,ISE_10,ISE_11,ISE_12,ISE_13,ISE_14,ISE_15))
mas<-t(cbind(mas_1,mas_2,mas_3,mas_4,mas_5,mas_6,mas_7,mas_8,mas_9,mas_10,mas_11,mas_12,mas_13,mas_14,mas_15))


for(b in 1:B){
    
    
    for(i in 1:15){
    x<-seq(min(mas[i,]),max(mas[i,]),by=0.01)
    x0<-sample(x,length(mas[i,]),replace=T)
    
    cross_validation_bandwidth<-optimize(LCV,interval=c(0,max(dist(mas))),mas=mas[i,])
    normal_scale_bandwidth<-1.06*min(sqrt(var(mas[i,])),IQR(mas[i,])/(qnorm(0.75)-qnorm(0.25)))*length(mas[i,])^(-1/5)
    density_CV<-normalkernel(x0,mas,cross_validation_bandwidth$minimum)
    density_normal<-normalkernel(x0,mas,normal_scale_bandwidth)
    density_theorical=densities[i,]
    for(m in 1:M){
    aux1[m]<-sum((density_CV-density_theorical)^2)
    aux2[m]<-sum((density_normal-density_theorical)^2)
    }
    ISE[b,i]<- mean(aux1)
    ISE[b,15+i]<-mean(aux2)
    }
}
    
