Thanks to Lindsey for supplying her code for HW2.


### S-PLUS CODE FOR A VARIABLE KERNEL DENSITY ESTIMATE 

## INPUT DATA
#data=input data 
data<-geyser$duration

## CHOICE OF BANDWIDTH
#Using Scale h
n<-length(data)
b<-(8/6)^{1/5}*sqrt(var(data))*n^{-1/5}

#Standardise the Kernel in Splus (sigma=0.37065, in Splus)
bopt<-b/0.37065

## CALCULATE AN ESTIMATE OF f
#bandwidth chosen by scale choice of h
#f=density estimate
f1<-ksmooth(data,kernel="normal",bandwidth=bopt)
f<-f1$y

## PLOT OF DATA 
b<-seq(0,6,0.5)
hist(data,breaks=b,ylim=c(0,0.6),prob=T,xlab="Duration",ylab="Density",main="P
robability Histogram of Data")

## PLOT OF f (WITH RUG PLOT INCLUDED)
plot(f1,type="l",xlab="Duration",ylab="Density",main="Pilot Estimate of f, Nor
mal Kernel")
rug(data)

## LEAST SQUARES CV - MINIMISE VIA GRID SEARCH
#f contains the pilot estimates 
#       matfi=matrix of f[Xi]
#       matfj=matrix of f[Xj]
fdat<-length(f)
matfi<-matrix(f,fdat,fdat)
matfj<-t(matfi)

#data contains the data points 
#       construct a matrix of data points (Xi-Xj)
#       mat[i,j]=abs(data[i]-data[j])

ndat<-length(data)
mat1<-matrix(data,ndat,ndat)
mat2<-t(mat1)
mat<-abs(mat1-mat2)

#construct a grid of values for h - vary this to find the minimum
hgrid<-seq(0.01,2.5,0.1)
ngrid<-length(hgrid)


# first term used in the calculation (K*K):
t1<-function(h,mat,matfi,matfj)
{
var<-(h*(1/sqrt(matfi)+1/sqrt(matfj)))^2  ##var of K*K
varv<-as.vector(var)                    ##convert to vector form to read
matv<-as.vector(mat)                    ##into dnorm
sum(dnorm(matv,mean=0,sd=sqrt(varv)))
}

#second term used in the calculation

t2<-function(h,matfj,mat)
{
math<-(sqrt(matfj)*mat)/h  
mathv<-as.vector(math)                ##convert to matrix form 
K<-dnorm(mathv,mean=0,sd=1)
Km<-matrix(K,ndat,ndat)              ##convert back to matrix form
mx<-((sqrt(matfj))/h)*Km
sum(mx)-sum(diag(mx))
}

## LEAST SQUARES CV - MINIMISE THIS VIA GRID SEARCH
Fh<-rep(0,ngrid)
for(i in 1:ngrid)
{
h<-hgrid[i]
Fh[i]<-((ndat^{-2})*t1(h,mat,matfi,matfj))-((2/ndat)*(1/(ndat-1))*t2(h,matfj,m
at))
}
# Plots Fh against values of h to find minimium
plot(hgrid,Fh,xlab="h",main="CV Curve for Geyser Data")

************************************************************************
## Generate standard normal variables then use the commands above 
## with various hgrids

data<-rnorm(100)
data<-sort(data)

## PLOT OF DATA 
b<-seq(-3,3,0.5)
hist(data,breaks=b,ylim=c(0,0.45),prob=T,xlab="x",ylab="Density",main="Probability Histogram of N(0,1) Data")

# Try this hgrid 
hgrid<-seq(0.1,2,0.1)   #plot shows min is between 0.3 and 0.5
hgrid<-seq(0.3,0.5,0.01)  #plot shows min is 0.41 

#Plot Fh against h to find minimum
plot(hgrid,Fh,xlab="h",main="CV Curve for N(0,1) Data")

## FINAL DENSITY ESTIMATE
#hstar=h* found from CV
hstar<-0.41

#term used in the final density estimate
calcf<-function(d)
{
x<-data-d
math<-(sqrt(f)*x)/hstar  
K<-dnorm(math,mean=0,sd=1)
mx<-sqrt(f)*K
sum(mx)
}

#Calculation of final density estimate
#DE=final density estimate using h*=0.27
DE<-rep(0,length(data))
for(i in 1:length(data))
{
d<-data[i]
DE[i]<-(1/(ndat*hstar))*calcf(d)
}

#plot final density estimate
plot(DE,xlab="",ylab="density",main="Final Density Estimate of f using variable bandwidth, h*=0.41")