{

#Let us suppose that the data is in xdata
x1=min(xdata)
x2=max(xdata)
delta=(x2-x1)/97
x1=x1-delta
x2=x2+delta
xeval=seq(x1,x2,delta)
nevals=length(xeval)

#xeval is the array of 100 points at which you want to estimate
#the PDFs

zz=sm.density(xdata,model="Normal",eval.points=xeval)
xdensityorig=zz$estimate	#PDF of the original data


xsimpdf=matrix(0,nrow=nsim,ncol=nevals)

#generate bootstrap samples of the data and estimate the PDF
#of each bootstrap sample..

for(i in 1:nsim){
x=sample(xdata,replace=T)
zz=sm.density(x,eval.points=xeval,model="Normal")
xsimpdf[i,]=zz$estimate
}

# the following will producethe boxplot of the PDFs along
# with the PDF of the original data..

xs=1:nevals
zz=boxplot(split(t(xevals),xs),plot=F,cex=1.0)
zz$names=rep("",length(zz$names))
z1=bxp(zz,ylim=range(xsimpdf,xdensityorig),xlab="",ylab="",style.bxp="old",cex=1.25)

z2=1:6
n1=1:6
z2[1]=z1[1]
z2[2]=z1[20]
z2[3]=z1[40]
z2[4]=z1[60]
z2[5]=z1[80]
z2[6]=z1[100]

n1[1]=xeval[1]
n1[2]=xeval[20]
n1[3]=xeval[40]
n1[4]=xeval[60]
n1[5]=xeval[80]
n1[6]=xeval[100]
n1=round(n1,dig=0)
n1=as.character(n1)

axis(1,at=z2,labels=n1,cex=1.00)
lines(z1,xdensityorig,lty=2,lwd=2)


#title(main="PDFs from Proportion Disagg. method")
}