### R commands for section 10.2.2 -- Choice of kernel
### G. Givens & J. Hoeting, Computational Statistics, Wiley, 2005.

## Normal kernel
f.normal <- function(x) { dnorm(x); }

## Uniform kernel
f.uniform <- function(x) { ifelse(abs(x)<1, 0.5, 0); }

## Epanechnikov kernel
f.epanechnikov <- function(x) { ifelse(abs(x)<1, (1-x^2)*3/4, 0); }

## Triangle kernel
f.triangle <- function(x) { ifelse(abs(x)<1, 1-abs(x), 0); }

## Biweight kernel
f.biweight <- function(x) { ifelse(abs(x)<1, (1-x^2)^2*15/16, 0); }

## Triweight kernel
f.triweight <- function(x) { ifelse(abs(x)<1, (1-x^2)^3*35/32, 0); }

## Compare different kernel functions
par(mfrow=c(2,2))
plot(c(-3,3),c(0,1.2),xlab="z", ylab="", main="Uniform Kernel", type="n")
curve(f.epanechnikov(x), c(-3,3), add=T, lwd=2, lty=1, col="black")
curve(f.normal(x), c(-3,3), add=T, lwd=2, lty=1, col="green")
curve(f.uniform(x), c(-3,3), add=T, lwd=2, lty=2, col="red")
legend(0.7,1,legend=c("Epanech","Normal"), lwd=c(2,2), lty=c(1,1), col=c("black","green"),cex=0.8)
plot(c(-3,3),c(0,1.2),xlab="z", ylab="", main="Triangle Kernel",type="n")
curve(f.epanechnikov(x), c(-3,3), add=T, lwd=2, lty=1, col="black")
curve(f.normal(x), c(-3,3), add=T, lwd=2, lty=1, col="green")
curve(f.triangle(x), c(-3,3), add=T, lwd=2, lty=2, col="red")
legend(0.7,1,legend=c("Epanech","Normal"), lwd=c(2,2), lty=c(1,1), col=c("black","green"),cex=0.8)
plot(c(-3,3),c(0,1.2),xlab="z", ylab="", main="Biweight Kernel",type="n")
curve(f.epanechnikov(x), c(-3,3), add=T, lwd=2, lty=1, col="black")
curve(f.normal(x), c(-3,3), add=T, lwd=2, lty=1, col="green")
curve(f.biweight(x), c(-3,3), add=T, lwd=2, lty=2, col="red")
legend(0.7,1,legend=c("Epanech","Normal"), lwd=c(2,2), lty=c(1,1), col=c("black","green"),cex=0.8)
plot(c(-3,3),c(0,1.2),xlab="z", ylab="", main="Triweight Kernel",type="n")
curve(f.epanechnikov(x), c(-3,3), add=T, lwd=2, lty=1, col="black")
curve(f.normal(x), c(-3,3), add=T, lwd=2, lty=1, col="green")
curve(f.triweight(x), c(-3,3), add=T, lwd=2, lty=2, col="red")
legend(0.7,1,legend=c("Epanech","Normal"), lwd=c(2,2), lty=c(1,1), col=c("black","green"),cex=0.8)

### Note: R function "density" provides choices of kernel functions as follows
### kernel = c("gaussian", "epanechnikov", "rectangular", "triangular", "biweight",
###            "cosine", "optcosine"),

## Example 10.5, Bimodal data, continued, page 293
## need data "bimodal.dat"
x = read.table("bimodal.dat")$V1

# Sheather-Jones bandwidth for normal kernel
bw.normal = bw.SJ(x)                   # Sheather-Jones bandwidth for normal kernel
#[1] 0.6881154
delta.normal = 1/(2*sqrt(pi))^(1/5)    # delta(K), Table 10.1 on page 292
bw.uniform = bw.normal*(9/2)^(1/5)/delta.normal     # equivalent bandwidth for uniform kernel, (10.32)
#[1] 1.19736
bw.epanechnikov = bw.normal*15^(1/5)/delta.normal   # equivalent bandwidth for Epanechnikov kernel
#[1] 1.523353
bw.triangle = bw.normal*24^(1/5)/delta.normal       # equivalent bandwidth for Triangle kernel
#[1] 1.673495
bw.biweight = bw.normal*35^(1/5)/delta.normal       # equivalent bandwidth for Biweight kernel
#[1] 1.804662
bw.triweight = bw.normal*(9450/143)^(1/5)/delta.normal      # equivalent bandwidth for Triweight kernel
#[1] 1.523353

## self-defined function to add "uniform" and "triweight" into "density" as "kernel"
density.self <- function(x, bw, kernel, n, from, to) {
  # variables follow R function "density"
  K <- function(x) { ifelse(abs(x)<1, 0.5, 0); };
  if(kernel=="triweight") {
    K <- function(x) { ifelse(abs(x)<1, (1-x^2)^3*35/32, 0); }
  };
  if(kernel=="epanechnikov") {
    K <- function(x) { ifelse(abs(x)<1, (1-x^2)*3/4, 0); }
  };
  if(kernel=="triangle") {
    K <- function(x) { ifelse(abs(x)<1, 1-abs(x), 0); }
  };
  if(kernel=="biweight") {
    K <- function(x) { ifelse(abs(x)<1, (1-x^2)^2*15/16, 0); }
  };
  if(kernel=="gaussian") {
    K <- function(x) { dnorm(x); }
  };
  xn = seq(from, to, length=n);        # x at which f(x) is calculated
  N = length(x);                       # number of observations
  temp = xn %*% t(rep(1,N)) - rep(1,n) %*% t(x);
  temp = K(temp/bw);
  y = apply(temp, 1, sum)/(N*bw);
  list(x=xn, y=y);
}

## plot Fig. 10.6 on page 294
par(mfrow=c(2,3))
temp = density.self(x, bw = bw.uniform, kernel="uniform", n=500, from=-2, to=18)
plot(c(-2,18),c(0,0.20),ann=F,frame.plot=T,type="n")
points(temp$x, temp$y, type="l", lty=1)
text(14, 0.17, "Uniform")
temp = density.self(x, bw = bw.epanechnikov, kernel="epanechnikov", n=500, from=-2, to=18)
plot(c(-2,18),c(0,0.20),ann=F,frame.plot=T,type="n")
points(temp$x, temp$y, type="l", lty=1)
text(14, 0.15, "Epanechnikov")
temp = density.self(x, bw = bw.triangle, kernel="triangular", n=500, from=-2, to=18)
plot(c(-2,18),c(0,0.20),ann=F,frame.plot=T,type="n")
points(temp$x, temp$y, type="l", lty=1)
text(14, 0.15, "Triangle")
temp = density.self(x, bw = bw.normal, kernel="gaussian", n=500, from=-2, to=18)
plot(c(-2,18),c(0,0.20),ann=F,frame.plot=T,type="n")
points(temp$x, temp$y, type="l", lty=1)
text(14, 0.15, "Normal")
temp = density.self(x, bw = bw.biweight, kernel="biweight", n=500, from=-2, to=18)
plot(c(-2,18),c(0,0.20),ann=F,frame.plot=T,type="n")
points(temp$x, temp$y, type="l", lty=1)
text(14, 0.15, "Biweight")
temp = density.self(x, bw = bw.triweight, kernel="triweight", n=500, from=-2, to=18)
plot(c(-2,18),c(0,0.20),ann=F,frame.plot=T,type="n")
points(temp$x, temp$y, type="l", lty=1)
text(14, 0.15, "Triweight")