# Density Estimation code
# This is the Old Faithful Eruption Duration data from the text

table1011 <- c(1.7, 1.9, 3.4, 3.8, 4.1, 4.5, 
1.7, 1.9, 3.5, 3.9, 4.1, 4.5, 1.7, 2.0, 3.5, 3.9, 4.1, 4.5,
1.7, 2.0, 3.5, 3.9, 4.1, 4.6, 1.7, 2.0, 3.5, 3.9, 4.2, 4.6,
1.8, 2.0, 3.6, 4.0, 4.2, 4.6, 1.8, 2.0, 3.6, 4.0, 4.3, 4.6,
1.8, 2.3, 3.7, 4.0, 4.3, 4.6, 1.8, 2.3, 3.7, 4.0, 4.3, 4.6,
1.8, 2.3, 3.7, 4.0, 4.3, 4.6, 1.8, 2.3, 3.7, 4.0, 4.3, 4.6,
1.8, 2.5, 3.7, 4.0, 4.3, 4.6, 1.8, 2.9, 3.7, 4.0, 4.4, 4.7,
1.8, 2.9, 3.7, 4.1, 4.4, 4.7, 1.9, 3.1, 3.8, 4.1, 4.4, 4.7,
1.9, 3.2, 3.8, 4.1, 4.4, 4.8, 1.9, 3.3, 3.8, 4.1, 4.4, 4.9,
1.9, 3.4, 3.8, 4.1, 4.5)

 # How to find the bandwidth being used in kernel density estimation

tab1011dens <- density(table1011)

> names(tab1011dens)
[1] "x"         "y"         "bw"        "n"         "call"      "data.name" "has.na"   
> tab1011dens$bw
[1] 0.3678302


 # How to use the bandwidth rule mentioned in the text

tab1011densNRD <- density(table1011,bw="nrd")
> tab1011densNRD$bw
[1] 0.4332222


 # These examples show how to change the default bandwidth rule 
 #  with the "adjust" argument 

hist(table1011,prob=T)
lines(density(table1011))

hist(table1011,prob=T,ylim=c(0,.6))
lines(density(table1011,adjust=.5))

hist(table1011,prob=T,ylim=c(0,.65))
lines(density(table1011,adjust=.25))


 # R code for the Figures

par(mfrow=c(2,1))
hist(table1011,prob=T,ylim=c(0,2.0),xlab="Old Faithful Eruption Duration",
main = quote("Rectangular kernel function, very small bandwidth " (Delta==0.01)))
lines(density(table1011,kernel="rectangular",bw=.01),col=4)
hist(table1011,prob=T,ylim=c(0,2.0),xlab="Old Faithful Eruption Duration",
main = quote("Rectangular kernel function, very large bandwidth "(Delta==1.5)) )
lines(density(table1011,kernel="rectangular",bw=1.5),col=4)
par(mfrow=c(1,1))


par(mfrow=c(3,1))
hist(table1011,prob=T,ylim=c(0,.60),xlab="Old Faithful Eruption Duration",
main = quote("Gaussian kernel function, bandwidth rule mentioned in the text " (Delta==0.433)))
lines(density(table1011,kernel="gaussian",bw=.433),col=4)
hist(table1011,prob=T,ylim=c(0,.60),xlab="Old Faithful Eruption Duration",
main = quote("Gaussian kernel function, automatic bandwidth rule " (Delta==0.368)))
lines(density(table1011,kernel="gaussian"),col=4)
hist(table1011,prob=T,ylim=c(0,.60),xlab="Old Faithful Eruption Duration",
main = quote("Gaussian kernel function, Sheather-Jones bandwidth rule " (Delta==0.185)))
lines(density(table1011,kernel="gaussian",bw="SJ"),col=4)
par(mfrow=c(1,1))