{

#simulate from a PDF, compute the stats for each Monte Carlo simulation
# and boxplot them along with those of the historical data

#Nonparametric Monte Carlo

data=matrix(scan("aismr.txt"),ncol=13,byrow=T)
data=matrix(scan("Leesferry-mon-data.txt"),ncol=13,byrow=T)
xdata=data[,6]

#suppose the historical data is in the vector 'xdata'
N=length(xdata)		#number of data points.


nsim=150	#number of simulations, each of length N

#points at which to estimate the PDF from the simulations and the observed

xeval=seq(min(xdata)-sd(xdata),max(xdata)+sd(xdata),length=100)
neval=length(xeval)
xsimpdf=matrix(0,nrow=nsim,ncol=neval)


meansim=1:nsim
mediansim=1:nsim
sdsim=1:nsim
skewsim=1:nsim
iqrsim=1:nsim

#generate Monte Carlo samples

library(sm)
band=hnorm(xdata)

for(i in 1:nsim){

#Nonparametric

#Bootstrap..
# you can replace this for monte carlo. E.g. x = rnorm(length(xdata), p1,p2)
# where p1 and p2 are the parameters of the normal distribution. Remember
# to comment out the boostrap command if using this..

xsim=sample(xdata,replace=T)

#Smooth Bootstrap from kernel density estimator...
#xsim=sample(xdata,replace=T)+band*rnorm(N)

#varkern is the variance of the kernel you choose..
#varkern=1		#variance of Normal kernel is 1
#denom=sqrt(1 + (band*band*varkern/var(xdata)))
#xsim=mean(xdata) + ((sample(xdata,replace=T) - mean(xdata) + band*rnorm(N))/denom)

xsimpdf[i,]= sm.density(xsim,eval.points=xeval,display="none")$estimate

meansim[i]=mean(xsim)
sdsim[i]=sd(xsim)
skewsim[i]=skew(xsim)
iqrsim[i]=diff(quantile(xsim,c(0.25,0.75)))
mediansim[i]=quantile(xsim,c(0.5))
}

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


#boxplot the stats. along with the the historical stats.

par(mfrow=c(1,5))

zz=boxplot(meansim,xlab="",ylab="Mean",plot=F)
zz$names=rep("",length(zz$names))
z1=bxp(zz,xlab="",ylab="Mean",cex=1.25)
points(z1,mean(xdata),col="red",cex=1.25)
title(main="Boxplot of Mean")

zz=boxplot(mediansim,xlab="",ylab="Mean",plot=F)
zz$names=rep("",length(zz$names))
z1=bxp(zz,xlab="",ylab="Mean",cex=1.25)
points(z1,quantile(xdata,c(0.5)),col="red",cex=1.25)
title(main="Boxplot of Median")

zz=boxplot(sdsim,xlab="",ylab="Mean",plot=F)
zz$names=rep("",length(zz$names))
z1=bxp(zz,xlab="",ylab="Mean",cex=1.25)
points(z1,sd(xdata),col="red",cex=1.25)
title(main="Boxplot of SD")

zz=boxplot(skewsim,xlab="",ylab="Mean",plot=F)
zz$names=rep("",length(zz$names))
z1=bxp(zz,xlab="",ylab="Mean",cex=1.25)
points(z1,skew(xdata),col="red",cex=1.25)
title(main="Boxplot of Skew")

zz=boxplot(iqrsim,xlab="",ylab="Mean",plot=F)
zz$names=rep("",length(zz$names))
z1=bxp(zz,xlab="",ylab="Mean",cex=1.25)
points(z1,diff(quantile(xdata,c(0.25,0.75))),col="red",cex=1.25)
title(main="Boxplot of IQR")

par(ask=TRUE)
par(mfrow=c(1,1))

xs=1:neval

#For R versions < 2.4.1
xsimpdf1=as.data.frame(t(xsimpdf))
zz=boxplot(split(xsimpdf1,xs),plot=F,cex=1.0)

#For R version 2.4.1
zz=boxplot(split(t(xsimpdf),xs),plot=F,cex=1.0)
zz$names=rep("",length(zz$names))
z1=bxp(zz,ylim=range(xsimpdf,xdensityorig),xlab="",ylab="",cex=1.25)

npts=10         #number of points to point..
n2=round(neval*(1:npts)/npts)

z2=z1[n2]
n1=xeval[n2]


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

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

title(main="PDFs from the simulations and the historical data")

}