Envelope Plots for fully specified Distribuions

Quantile Confidence Bands

required libraries

library(ggplot2)    # for graphics

Numerical calculation of quantile function

qF <-
function(p,F,np=500) {
      bisect<-function(x,low,high,acc=0.0001) {
            repeat {
                mid<-(low+high)/2
                if(F(mid)<x) low<-mid
                else high<-mid
                if(high-low<acc) break
            }
            mid
      }
      if(missing(p)) p<-(1:np)/(np+1)
      else {
          np<-length(p)
          p<-sort(p)
      }    
      y<-0*p
      minx<-0
      if(F(minx)>p[1]) {
            repeat {
                minx<-minx-1
                if(F(minx)<p[1]) break
            }
       }
       maxx<-0     
       repeat {
            maxx<-maxx+1
            if(F(maxx)>p[np]) break
       }      
       y[1]<-bisect(p[1],minx,maxx)
       if(np>1)
          for(i in 2:np) y[i]<-bisect(p[i],y[i-1],maxx)
      
      y
}

Envelope Plot

envelopePlot <-
function (x, pF,alpha=0.05,Out=FALSE) 
{
    if(length(x)<100) n<-length(x)
    else {
        n<-100
        x<-quantile(x,c(1:100)/101)
    }    
    y <- ppoints(n)
    if(alpha==0.1) alpha<-exp(-2.706-0.516*log(n))
    if(alpha==0.05) alpha<-exp(-3.37-0.544*log(n))
    if(alpha==0.01) alpha<-exp(-5.11-0.551*log(n))
    lower <- qbeta(alpha[1]/2, 1:n, n - 1:n + 1)
    upper <- qbeta(1 - alpha[1]/2, 1:n, n - 1:n + 1)
    qy<-qF(F=pF) 
    qFun<-function(p) {approx((1:500)/501,qy,xout=p)$y} 
    plt<- qplot(sample = x, geom = 'blank') + stat_qq(distribution = qFun,colour="red")        
    dtaLim<-data.frame(x=qFun(y),y=qFun(lower))
    dtaLim<-dtaLim[!is.na(dtaLim[,2]), ]
    plt<-plt+geom_line(aes(x,y),colour="blue",data=dtaLim,size=1.05)
    dtaLim<-data.frame(x=qFun(y),y=qFun(upper))
    dtaLim<-dtaLim[!is.na(dtaLim[,2]), ]
    plt<-plt+geom_line(aes(x,y),colour="blue",data=dtaLim,size=1.05)
    if(Out) return(plt)
    print(plt)
    
}

Example

x <- rexp(100)
envelopePlot(x, pexp)