simboxpdf = function(integerMonth,simflow,obsflow)
{

library(sm)

test = matrix(obsflow, ncol=12,byrow=T)

imon = integerMonth

data=test[,imon]
#data=log(data)
xx=data

#bandd=hsj(xx)		#sheather Jones approach for bandwidth
bandd=hnorm(xx)		#for the Normal bandwidth..
#bandd = bandd*2.214    #for EP kernel
#band=bandd

#go one bandwidth from the min and max of the data..
xlow=min(data)-bandd
#xlow=max(0,xlow)
xhigh=max(data)+bandd

#create 50 points equally spaced between xlow and xhigh..
xeval=seq(xlow,xhigh,length=50)
neval=length(xeval)

zz=sm.density(data,eval.points=xeval,display="none")

# When you estimate PDF in the log space you have to divide by
# exp(of evaluation points) to get the pdf in the original space
# See Lall et al. (1996) paper for the formula..

#pdf of the observed data..
#xpdfo=zz$estimate/exp(xeval)
xpdfo=zz$estimate

#now read the simulations..
xsim=matrix(scan(simflow), ncol=nsim, byrow=T)
		
xpdf=matrix(0,ncol=nsim, nrow=neval)

#get pdf of the simulations..

for(i in 1:nsim){
	xx=xsim[,i]
	#xx=log(xx)

	#bandd=hsj(xx)
	#bandd=hnorm(xx)
	#bandd = bandd*2.214     #for EP kernel
	#band=bandd

	zz=sm.density(xx,eval.points=xeval,display="none")
	#xpdf[,i]=zz$estimate/exp(xeval)
	xpdf[,i]=zz$estimate		
}


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

xs=1:neval

yrange=c(yrange,(range(xpdf,xpdfo)))
assign("yrange", yrange, env = .GlobalEnv)

par(col.lab="black",cex.lab=1.25,cex=1)

xpdf=as.data.frame(xpdf)
zz=myboxplot(split(xpdf,xs),plot=F,cex=.75)
zz$names=rep(" ",length(zz$names))
z1=bxp(zz,ylim=range(xpdf,xpdfo),xlab=paste("flow (cms)"),ylab="PDF",cex=.75,axes=F)
box(lty="solid",lwd=2)

z2=1:10
n1=1:10
n2=round(neval/10)
z2[1]=z1[1]
z2[2]=z1[2*n2]
z2[3]=z1[3*n2]
z2[4]=z1[4*n2]
z2[5]=z1[5*n2]
z2[6]=z1[6*n2]
z2[7]=z1[7*n2]
z2[8]=z1[8*n2]
z2[9]=z1[9*n2]
z2[10]=z1[neval]

#xeval=exp(xeval)

n1[1]=xeval[1]
n1[2]=xeval[2*n2]
n1[3]=xeval[3*n2]
n1[4]=xeval[4*n2]
n1[5]=xeval[5*n2]
n1[6]=xeval[6*n2]
n1[7]=xeval[7*n2]
n1[8]=xeval[8*n2]
n1[9]=xeval[9*n2]
n1[10]=xeval[neval]

n1=round(n1,dig=0)
n1=as.character(n1)

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

}