# For all these examples, nsim is way too small.
MEAN <- 0 # mean of zero
SIGMA <- 100 # standard deviation of 100
PAR <- vec2par(c(MEAN,SIGMA),type='nor') # make parameter object
F <- c(0.5, 0.8, 0.9, 0.96, 0.98, 0.99) # nonexceed probabilities
# nsim is small for speed of example not accuracy.
CI <- genci(PAR,n=10,F=F,nsim=20)
plot(CI$nonexceed_prob,CI$true,type='l',lwd=2)
lines(CI$nonexceed_prob,CI$lower,col=2)
lines(CI$nonexceed_prob,CI$upper,col=3)
pdf("twoCIplots.pdf")
# The qnorm call has been added to produce "normal probability"
# paper on the horizonal axis. The parent is heavy-tailed.
GEV <- vec2par(c(10000,1500,-.3),type='gev') # a GEV distribution
CI <- genci(GEV,n=20,nsim=200,edist='gno')
ymin <- log10(min(CI$lower[! is.na(CI$lower)]))
ymax <- log10(max(CI$upper[! is.na(CI$upper)]))
plot( qnorm(CI$nonexceed_prob),log10(CI$true),type='l',
ylim=c(ymin,ymax),lwd=2)
lines(qnorm(CI$nonexceed_prob),log10(CI$lower),col=2)
lines(qnorm(CI$nonexceed_prob),log10(CI$upper),col=3)
# another error distribution model
CI <- genci(GEV,n=20,nsim=200,edist='aep4')
lines(qnorm(CI$nonexceed_prob),log10(CI$lower),col=2, lty=2)
lines(qnorm(CI$nonexceed_prob),log10(CI$upper),col=3, lty=2)
dev.off()
Run the code above in your browser using DataLab