{

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

library(MASS)

data=matrix(scan("aismr.txt"),ncol=13,byrow=T)
xdata=data[,11]

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

#Monte Carlo...
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


for(i in 1:nsim){

#parametric..
#simulate from Lognormal
#replace with any other distribution..
xsim=rlnorm(N,meanlog=mean(log(xdata)), sdlog=sd(log(xdata)))

#Gamma..
#zparam=fitdistr(xdata,"gamma")
#xsim=rgamma(N,shape=zparam$estimate[1],scale=1/zparam$estimate[2])


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

#estimate the PDF at the evaluation points based on the simulated data
#replace with the appropriate distribution..
xsimpdf[i,]=dlnorm(xeval,meanlog=mean(log(xsim)), sdlog=sd(log(xsim)))

zz=fitdistr(xsim,"gamma")
xsimpdf[i,]=dgamma(xeval,shape=zz$estimate[1],scale=1/zz$estimate[2])

}

#PDF of the original data..
#replace with the appropriate distribution..
xdensityorig=dlnorm(xeval,meanlog=mean(log(xdata)), sdlog=sd(log(xdata)))

#xdensityorig=dgamma(xeval,shape=zparam$estimate[1],scale=1/zparam$estimate[2])

#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(log(xdata)),col="red",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)
#boxplot the PDFS
par(mfrow=c(1,1))
xs=1:neval


# For R version < 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 plot on the x-axis..
n2=round(neval*(1:npts)/npts)

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


n1=round(n1,dig=2)
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")


}