library(sm)

# pass in your data in the variable data. If not you can scan
#it in
#rain=matrix(scan("aismr.txt"),ncol=13,byrow=T)
#data=rain[,13]

#data from exponential distribution to test the boundary transformation..
data=rexp(250)		
data1=data

n=length(data)

#add the PDF from a normal kernel
zz=sm.density(data,display="none")
xeval=zz$eval.points	#points at which the PDF is evaluated. You can supply your own
xeval=seq(min(data), max(data)+zz$h, length=500)
				#if you like
band=zz$h	#reference bandwidth for the normal kernel

#band=hsj(data)	# to get Sheather-Jones bandwidth using Normal Kernel
#band=hcv(data)	# to get Least Squares Cross Validation Bandwidth
			# using Normal kernel# 
# to get estimates of pdf with Normal kernel with a different bandwidth
#zz=sm.density(data, hsj(data))
#zz=sm.density(data, hcv(data))
#zz=sm.density(data, hnorm(data))	# the default
# etc..


#bandwidth for Bisquare kernel..
band=band*2.623
neval=length(xeval)
pdfbs=1:neval

#kernel estimator using a Bisquare kernel
for(i in 1:neval){
xx=(xeval[i]-data)/band
xx[abs(xx) >= 1]=1		# if the t value is beyond -1,1 replace it
				# with 1 (i.e. making the weights 0.
#For normal kernel the above line should be commented..

yy=(15./16.) * (1. - xx*xx)^2		#if you want to use EP kernel
				# you replace this with yy=(3/4) * (1. - xx*xx)
pdfbs[i]=sum(yy)/(n*band)
}


#If you require boundary correction..

#get the bandwidth of the log data
data=log(data)
zz=sm.density(data,display="none")
band=zz$h
#band=hsj((data))		#if you want the Sheather Jones bandwidth

#xevallog=zz$eval.points
#xevallog=xevallog[xevallog > 0]	#estimate only at points > 0.

pdfbsbc=1:neval

#bandwidth for Bisquare kernel..
band=band*2.623

#estimate the PDF from the Bisquare kernel
for(i in 1:neval){
xx=(log(xeval[i])-data)/band
xx[abs(xx) >= 1]=1		# if the t value is beyond -1,1 replace it
yy=(15./16.) * (1. - xx*xx)^2
#yy[abs(xx) >= 1]=0
pdfbsbc[i]=sum(yy)/(n*band*xeval[i])
}



## now plot..
hist(data1,probability=T,xlab="Rain", ylab="PDF",ylim=range(pdfbs,pdfbsbc),xlim=range(xeval))

#boundary corrected PDF from bisquare..
lines(xeval, pdfbsbc, col="red")

#no boundary correction - PDF from bisquare kernel.
lines(xeval, pdfbs, col="blue")

#exponential density fitted to the data..
lines(xeval,dexp(xeval,1/mean(data1)), col="green")

yy=rep(min(pdfbs,pdfbsbc),n)
points(data1,yy,col='red')
#normal kernel - no boundary correction..
sm.density(data, add=T, lwd=2)