## Slide 5
d <- density(faithful$waiting, kernel="rectangular", bw=2/sqrt(3))
plot(d, xlab="Waiting times (min)", main="h=2")

## Gaussian figure (slide 8)
x <- c(1, 3, 5, 6, 6.2, 6.6, 7, 8, 10)
n <- length(x)
xx <- seq(-2, 10, len=101)
hist(x, border="white", main="", freq=FALSE)
rug(x, col="red", lwd=3)
for (i in 1:n) lines(xx, dnorm(xx, x[i], 0.75)/10, col="#008DFFB2", lwd=3)
d <- density(x, bw=0.75)
lines(d, lwd=3)

## Kernels (slide 10)
n <- 300
xx <- seq(-3,3,len=n)
df <- rbind(data.frame(x=xx, Density=dnorm(xx), kernel="Gaussian"),
            data.frame(x=xx, Density=3/4*(1-xx^2)*(abs(xx) <= 1), kernel="Epanechnikov"),
            data.frame(x=xx, Density=(1 - abs(xx)^3)^3 * (abs(xx) <= 1), kernel="Tri-cube"))
xyplot(Density~x,df,group=kernel,auto.key=list(points=FALSE,lines=TRUE,columns=3),type="l")

## Changing bw (slide 11)
x <- faithful$waiting
par(mfrow=c(1,3))
plot(density(x, bw=10), main="h=10", xlab="Waiting time (min)", lwd=3)
plot(density(x, bw=2.66), main="h=2.66", xlab="Waiting time (min)", lwd=3)
plot(density(x, bw=1), main="h=1", xlab="Waiting time (min)", lwd=3)
bw.ucv(x)

## Different bandwidth rules
par(mfrow=c(1,2))
d <- density(x, bw="nrd")
plot(d, main=paste("Normal: h = ", round(d$bw,2), sep=""), xlab="Waiting time (min)", lwd=3)
d <- density(x, bw="ucv")
plot(d, main=paste("CV: h = ", round(d$bw,2), sep=""), xlab="Waiting time (min)", lwd=3)

## Different kernels
par(mfrow=c(1,3))
plot(density(x, kernel="gaussian", bw="ucv"), main="Gaussian", xlab="Waiting time (min)", lwd=3)
plot(density(x, kernel="epanechnikov", bw="ucv"), main="Epanechnikov", xlab="Waiting time (min)", lwd=3)
plot(density(x, kernel="biweight", bw="ucv"), main="Biweight", xlab="Waiting time (min)", lwd=3)

## 2D (slides 23-25)
require(KernSmooth)
par(mfrow=c(1,1))
with(faithful, plot(eruptions,waiting,pch=19,xlab="Duration (min)",ylab="Waiting (min)",col="blue"))
fit <- bkde2D(faithful, bandwidth=c(0.7, 7))
contour(fit$x1, fit$x2, fit$fhat,xlab="Duration (min)",ylab="Waiting (min)")
myPal <- colorRampPalette(c("white", "blue"))
filled.contour(fit$x1, fit$x2, fit$fhat,xlab="Duration (min)",ylab="Waiting (min)", color.palette=myPal)
persp(fit$x1, fit$x2, fit$fhat, shade=0.15, col="#BFE6FF", theta=30, phi=25, d=5, ticktype="detailed", xlab="Duration (min)",ylab="Waiting (min)",zlab="Density")
require(rgl)
persp3d(fit$x1, fit$x2, fit$fhat, col="#BFE6FF", xlab="Duration (min)",ylab="Waiting (min)",zlab="Density")

## Bootstrap (slides 29-30)
## Illustration
require(boot)
xx <- seq(32,107,len=101)
fun <- function(x, ind, xx) {
  d <- density(x[ind])
  approx(d$x, d$y, xx)$y
}
out <- boot(faithful$waiting, fun, R=50, xx=xx)
d <- density(faithful$waiting)
plot(d,main="",xlab="Waiting time (min)",ylim=c(0,.041))
matplot(xx,t(out$t),add=TRUE,lty=1,col="gray",type="l")
lines(d,col="red",lwd=3)

## Pointwise band
xx <- seq(35, 100, len=101)
out <- boot(faithful$waiting, fun, R=999, xx=xx)
l <- apply(out$t, 2, quantile, p=.025)
u <- apply(out$t, 2, quantile, p=.975)

## Global band
a1 <- exp(seq(log(.01), log(.00001), len=49))
a2 <- numeric(length(a1))
fb <- out$t
for (i in 1:length(a1)) {
  lower <- apply(fb, 2, quantile,p=a1[i]/2)
  upper <- apply(fb, 2, quantile,p=1-a1[i]/2)
  a2[i] <- mean(apply(t(fb) < lower,2,any) | apply(t(fb) > upper,2,any))
}
plot(a1,a2,xlim=rev(range(a1)),type="l",xlab=expression(alpha[pw]),ylab=expression(alpha))
abline(h=.05,col="red")
alpha <- max(a1[a2 < .05])
L <- apply(fb, 2, quantile, p=alpha/2)
U <- apply(fb, 2, quantile, p=1-alpha/2)

## Plot w/ global, local bands
plot(d, main="", xlab="Waiting time (min)", ylim=c(0,.045))
polygon(c(xx,rev(xx)), c(L,rev(U)), col="gray90", border=FALSE)
polygon(c(xx,rev(xx)), c(l,rev(u)), col="gray80", border=FALSE)
lines(d, lwd=2.5)