#-------------------------------------------------------------------
calcdens <- function(aphat,dphat) {

#  This function computes the density of p at the desired p
#  points (dphat).  The actual points (aphat) are used to
#  compute the density.


     #  compute the optimal bandwidth for the density estimation
     #  (see Silverman for this formula)
  
     R <- quantile(aphat,c(.75)) - quantile(aphat,c(.25))
     hd1.old <- 0.79 * R * (length(aphat)**(-0.2))

     k2 <- 0.1428572
     k3 <- 0.7142857
     hd1 <- (k2**(-2/5))*(k3**(1/5))*(0.212**(-1/5))* (R/1.34)*(length(aphat)**(-1/5)) 
    # print(hd1)
     # loop through data to get the kernel density estimate at each point


     fhat <- vector("numeric",length(dphat))
     fhat.old <- vector("numeric",length(dphat))

     for (i in 1:length(dphat)) {
	  xo <- dphat[i]
	  arg <- (aphat-xo)/hd1 
	  arg2 <- (aphat-xo)/hd1.old 
	  kval <- (15/16)*((1-(arg**2))**2)*(abs(arg)<=1)
	  kval2 <- dnorm(arg2,0,1)*(arg2<=1)
	  fhat[i] <- (1/(length(aphat)*hd1))* sum(kval)
	  fhat.old[i] <- (1/(length(aphat)*hd1.old)*sum(kval2))
     }
	 
     dens1 <- fhat
     dens1
}