######################################################################################
###################### KERNEL DENSITY ESTIMATION #####################################
######################################################################################


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


#-- KERNEL ESTIMATOR

kern_est<-function(x,h,mas){

# Description:
#   This function implements the kernel estimator expression:
#       f.hat(x)=1/(n*h) sum_{i=1}^{n}{K((x-X_i))/h},
#   for a Normal kernel.

# Arguments:
#   x<- The point(s) where to evaluate the estimator.
#   h<- The bandwidth to use.
#   mas<- The sample X_1,...,X_n to use.    

    return(1/(length(mas)*h)*apply(dnorm(outer(x,mas,"-")/h),1,sum))
}


#-- LCV

LCV<-function(h,mas,tri_matrix){

# Description:
#   This function implements the the least cross validation expression:

#       LCV(h)= 1/(n^{2} * h)*sum_{i,j}{K*K ((X_i-X_j)/h)} 
#             - 4/(n*(n-1)*h) *sum_{i=1}^{n} * sum_{j>i}{K((X_i-X_j)/h)},

#   for a Normal kernel (so the convolution K*K corresponds to a Normal density whith mu=0 and sd=sqrt(2)).
#   When minimazing this function over h, the cross validation bandwidth is obtained.

# Arguments:
#   h<- The bandwidth to use.
#   mas <- The sample X_1,...,X_n to use.
#   tri_matrix<- An upper triangular matrix (logical) of dimension (n x n).

    return(sum(dnorm(outer(mas,mas,"-")/h,sd=sqrt(2)))/(length(mas)^2*h)
    - (4/(length(mas)*(length(mas)-1)*h)) * sum(dnorm((outer(mas,mas,"-")/h)[tri_matrix])))
}


#-- ISE INTEGRAND

ISE_int<-function(ap_dens,th_dens){

# Description:
#   This function returns the expression:
#       (f.hat(x)-f(x))^2,
#   where f.hat() is the kernel estimator and f() is the theorical density.

# Arguments:
#   ap_density<- The approximated density (kernel estimator evaluated over a(-) certain point (s).)
#   th_density<- The theorical density (evaluated a(-) certain point (s).)

    return((ap_dens-th_dens)^2)
}


#-- TRAPEZIUM RULE

trap_rule<-function(h,th_dens,est){

# Description:
#   This functions returns the integral aproximation of 

#           Int(g)=ISE(h)=int_{a}^{b}((f.hat(x;h)-f(x))^2 dx),

#   where f.hat is the kernel estimator with bandwidth h of the true density f,
#   using the trapezium rule. This rule consists on aproximating I(g) by

#           h/2 * (g(a)+g(b)) + h* sum(g(grid[2:n-1])),

#   where (a=grid[1],grid[2],...,grid[n]=b) is a partition of the interval [a,b].
#   Note that for this case, the function g() is ISE(), which has as arguments 
#   f.hat and f. 

# Arguments:
#   h: The bandwidth h.
#   th_dens: The values of the real function evaluated over a certain grid.
#   est: The values of the estimated function over a certain grid (the same as the previous one).


# Output:
#   This function returns an approximation of I(g) (a value).

            
    return(h/2 * (ISE_int(est[1],th_dens[1]) + ISE_int(est[length(est)],th_dens[length(th_dens)])) +
    h * sum(ISE_int(est[2:(length(est)-1)],th_dens[2:(length(th_dens)-1)])) )
}


#############################################################################################
################################### MAIN FUNCTION ###########################################
#############################################################################################


main<-function(M,grid,grid_opt,B,n,tri_matrix){

# Description:
#   This function computes both cross-validation and normal-scale bandwiths for a sample of 
#     a certain density. In addition, the kernel estimator is performed whith this bandwidths,
#   and the ISE criterion is computed.

# Arguments:
#   M<- A list including the sample generating function, its corresponding density function and
#       the name of the density.
#   grid<- A grid where the trapezium rule is carried out.
#   grid_opt<- Subintervals where the optimization of LCV is carried out.
#   B<- Number of simulations.
#   n<- Sample length.
#   tri_matrix<- An upper triangular matrix (logical) of dimension (n x n).


# Output:
#   This function returns the bandwidths (h) and the corresponding ISE (ISE) for B simulations.
#   It also returns the first five graphics associated to the LCV function (over a certain
#     range of bandwidths) and a comparative graphic for both theorical and approximated densities (for
#   both cross validation and normal scale bandwidths) over a grid. In addition, this function returns
#   the name of the Marron and WandŽ density that is being computed (M[[3]]).

    set.seed(0)

    for(b in 1:B){
            
        # 1. Extract the sample generating function from M and obtain the sample.

        mas<-M[[1]](n)
        
        # 2. Compute the cross - validation bandwidth by minimazind the LCV function (in subintervals).
        #    Plot the least cross validation function (for the first five simulations). 
        
        #if(b<=5){
        #   for (i in 1:length(grid[grid>0])){CV[i]= LCV(h=grid[grid>0][i],mas)}
        #   plot(grid[grid>0],CV,type="l",xlab="h",ylab="LCV(h)",main=M[[3]])
        #   windows()
        #}  

        opt<-apply(grid_opt,1,function(x) as.numeric(optimize(LCV,interval=c(x[1],x[2]),mas=mas,tri_matrix=tri_matrix)))

        h[b,1]<-opt[1,][which(opt[2,]==min(opt[2,]))]                           # h_cv
    
        # 3. Compute the normal - scale bandwidth
    
        h[b,2]<-1.06*min(sd(mas),IQR(mas)/(qnorm(0.75)-qnorm(0.25)))*length(mas)^(-1/5) # h_normal

        # 4. Extract the the density function from matrix_mas_density       

        th_dens<-M[[2]]
        
        # 5. Compute the associated ISE. Plot a comparative graphic for both theorical and approximated density
        #    (for the first five simulations)

        est_cv<-kern_est(grid,h[b,1],mas)       # Kernel estimation over grid with h_cv
        est_normal<-kern_est(grid,h[b,2],mas)   # Kernel estimation over grid with h_normal
        th_dens<-th_dens(grid)              # Theorical density evaluated over grid
        
        #if(b<=5){
        #   plot(th_dens~grid,type="l",xlab="",ylab="",main=M[[3]],sub="Cross Validation Bandwidth",ylim=c(0,max(th_dens)+0.5))
        #   lines(est_cv~grid,lty=3)
        #   windows()
        #   plot(th_dens~grid,type="l",xlab="",ylab="",main=M[[3]],sub="Normal Scale Bandwidth",ylim=c(0,max(th_dens)+0.5))
        #   lines(est_normal~grid,lty=3)
        #   windows()

        #}


        ISE[b,1] <- trap_rule(h[b,1],th_dens,est_cv)        # ISE(h_cv) 
        ISE[b,2] <- trap_rule(h[b,2],th_dens,est_normal)    # ISE(h_normal)

        
    }

    return(list(h,ISE,M[[3]]))

}


###############################################################################################
#################################### SIMULATION ###############################################
###############################################################################################


library(nor1mix)

#-- List of the sample generating functions, its corresponding theorical densities and its names (Marron
#   and Wand, 1992).

M_1<-c(function(x) rnorMix(x,MW.nm1),function(x)dnorMix(x,MW.nm1),"Gaussian density")
M_2<-c(function(x) rnorMix(x,MW.nm2),function(x)dnorMix(x,MW.nm2),"Skewed Unimodal Density")
M_3<-c(function(x) rnorMix(x,MW.nm3),function(x)dnorMix(x,MW.nm3),"Strongly Skewed Density")
M_4<-c(function(x) rnorMix(x,MW.nm4),function(x)dnorMix(x,MW.nm4),"Kurtotic Unimodal Density")
M_5<-c(function(x) rnorMix(x,MW.nm5),function(x)dnorMix(x,MW.nm5),"Outlier Density")
M_6<-c(function(x) rnorMix(x,MW.nm6),function(x)dnorMix(x,MW.nm6),"Bimodal Density")
M_7<-c(function(x) rnorMix(x,MW.nm7),function(x)dnorMix(x,MW.nm7),"Separated Bimodal Density")
M_8<-c(function(x) rnorMix(x,MW.nm8),function(x)dnorMix(x,MW.nm8),"Asymmetric Bimodal Density")
M_9<-c(function(x) rnorMix(x,MW.nm9),function(x)dnorMix(x,MW.nm9),"Trimodal Density")
M_10<-c(function(x) rnorMix(x,MW.nm10),function(x)dnorMix(x,MW.nm10),"Claw Density")
M_11<-c(function(x) rnorMix(x,MW.nm11),function(x)dnorMix(x,MW.nm11),"Double Claw Density")
M_12<-c(function(x) rnorMix(x,MW.nm12),function(x)dnorMix(x,MW.nm12),"Asymmetric Claw Density")
M_13<-c(function(x) rnorMix(x,MW.nm13),function(x)dnorMix(x,MW.nm13),"Asym. Db. Claw Density")
M_14<-c(function(x) rnorMix(x,MW.nm14),function(x)dnorMix(x,MW.nm14),"Smooth Comb Density")
M_15<-c(function(x) rnorMix(x,MW.nm15),function(x)dnorMix(x,MW.nm15),"Discrete Comb Density")

M<-list(M_1,M_2,M_3,M_4,M_5,M_6,M_7,M_8,M_9,M_10,M_11,M_12,M_13,M_14,M_15)


#-- Number of simulations

B<-500

#-- Grid (for the trapezium rule and also for the LCV graphic)

grid<-seq(-4,4,by=0.05)

#-- Subintervals (for the LCV optimization)

grid_opt<-cbind(c(0,seq(0,4,by=0.5)),c(seq(0,4,by=0.5),0))[-1,][-length(seq(0,4,by=0.5)),]

#-- Length of the samples.

n=100

#-- Inicializations

CV<-numeric(length(grid[grid>0]))           # Vector where to place the LCV evaluations.
tri_matrix<-upper.tri(matrix(0,nr=n,nc=n))  # Upper triangular matrix (logical) of dimension (n x n). 
h<-matrix(0,ncol=2,nrow=B)                  # Matrix of bandwidths (first column: h_cv; second column: h_normal) 
ISE<-matrix(0,ncol=2,nrow=B)                # Matrix of ISEs (first column: ISE(h_cv); second column: ISE(h_normal))


#-- Results of simulations

results<-lapply(M,main,grid=grid,grid_opt=grid_opt,B=B,n=n,tri_matrix=tri_matrix)


###############################################################################################
######################################### ANALYSIS ############################################
###############################################################################################


#-- Descriptive analysis of the results for each density.

minimums<-lapply(results,function(x) lapply(x[-3],function(x) apply(x,2,min)))
maximums<-lapply(results,function(x) lapply(x[-3],function(x) apply(x,2,max)))
means<-lapply(results,function(x) lapply(x[-3],function(x) apply(x,2,mean)))
medians<-lapply(results,function(x) lapply(x[-3],function(x) apply(x,2,median)))
sds<-lapply(results,function(x) lapply(x[-3],function(x) apply(x,2,sd)))

#-- Percentage of times ISE(h_normal)<ISE(h_cv) (i.e, estimation with h_normal is better than with h_cv)
#   for each density.

percentages<-lapply(results,function(x) sum(apply(x[[2]],1,function(x) x[1]>x[2]))/length(x))

#-- Boxplots of ISE(h_cv) and ISE(h_normal) (by pairs) for each density.

lapply(results,function(x) {boxplot(x[[2]],names=c("h_cv","h_ns"),main=x[[3]]);windows()})
