Qc<-function(x,quantile=0.5){ 
  #return sample quantile when ties are present based on the Mid CDF
          x<-sort(x)
          n<-length(x)
          qx<-umid(x)
          value<-NULL
         
          for(i in 1:(n-1))
           if(quantile<qx[n]){
            if(quantile==qx[i])
             value=x[i]
            if(quantile>qx[i]&&quantile<=qx[i+1]){
             b=(x[i+1]-x[i])/(qx[i+1]-qx[i])  #slope of regression
             a=x[i]-b*qx[i]                   #intercept of regression
             value=a+b*quantile}
             }
         if(quantile<=qx[1]){ 
                value=x[1]
                }
          if(quantile>=qx[n]){
                   value=x[n]
                }
          return(value)
        }



rd<-function(x,d=2) return(round(x,d))

P_<-function(p,n){  #converts quantiles(j-1)/(n-1) to (j-.5)/n
        val<-p+(p-.5)/(n-1)
        return(val)}

#Given P returns predict quantile (Normal)from Inverse 
#of mid-dist Function

qhat<-function(x,CI=.975){
        x<-sort(x)
        n<-length(x)
        p<-NULL
        for(k in 1:n)
                p[k]<-(k-.5)/n
        val<-NULL
        u<-NULL
        for(g in 1:n)
                u[g]<-p[g]+sqrt((p[g]*(1-p[g]))/n)*qnorm(CI)
        for(i in 1:n)
                val[i]=Qc(x,u[i])
        return(val)
}


Ftild<-function(x){        
        x<-sort(x)
        seg<-NULL        
        for(i in 1:(length(x)-1)){        
                count=1                
                for(j in i:(length(x)-1)){
                        if(i==1||x[i]>x[i-1]){
                                if(x[j]==x[j+1]){
                                        count=count+1        
                                }                                

                                else{
                                        break
                                }
                        }                        
                }                
                        if(i==1||x[i]>x[i-1])
                        seg=c(seg,count)
        }                
        if(x[length(x)]!=x[length(x)-1])
                        seg=c(seg,count)
        return(seg)
}

umid<-function(x){  #return the Mid CDF value
        x<-sort(x)
        rep<-Ftild(x)
        val<-rep[1]/2
        n<-length(rep)
        nx<-length(x)
        rt<-0
        for(i in 2:n){
                rt=rt+rep[i-1]
                val[i]=rt+rep[i]/2
                }
        val<-val/nx
        v2<-val[1]
        index<-2
        for(j in 2:nx)
                if(x[j]==x[j-1])
                        v2[j]=v2[j-1]
                else{                
                        v2[j]=val[index]
                        index=index+1
                }
        return(v2)
}

qnfplot<-function(x){       #Normal QQplot 
        x<-sort(x)
        u<-umid(x)
        n<-length(x)
        phi<-qnorm(P_(u,n))
        plot(phi,x,,xlab="theretical Normal quantiles",
                ylab="data quantiles",main="Normal QQplot with 95%CI",pch=20)
        points(qnorm(u),qhat(x,.975),col="blue",type="l")
        points(qnorm(u),qhat(x,.025),col="blue",type="l")
}


qufplot<-function(x){       #Uniform QQplot
        x<-sort(x)
        u<-umid(x)
        plot(u,x,,xlab="theretical Uniform quantiles with 95%CI",
		ylab="data quantiles",main="Uniform QQplot",pch=20)
        points(u,qhat(x,.975),col="blue",type="l")
        points(u,qhat(x,.025),col="blue",type="l")
}

qiq<-function(x){
  #returns vector standardized according to the MidQuartiles and the 2IQR
        q3<-Qc(x,.75)
        q1<-Qc(x,.25)
        mq<-(q3+q1)/2
        sq<-2*(q3-q1)
        xn<-(x-mq)/sq
        return(xn)}

qiq1<-function(x,p){  
   #returns a specific value that is standardized according to the 
   #MidQuartiles and the 2IQR
        q3<-Qc(x,.75)
        q1<-Qc(x,.25)
        mq<-(q3+q1)/2
        sq<-2*(q3-q1)
        xn<-(Qc(x,p)-mq)/sq
        return(xn)}

qiqn<-function(u,x){ #returns a vector of standardize normals (MD and 2IQR) 
                #with mean=mean(x) and sigma=stddev(x)
        n<-length(x)
        xb<-mean(x)
        sx<-sqrt(var(x))
        q1<-qnorm(.25,xb,sx)
        q3<-qnorm(.75,xb,sx)
        md<-(q3+q1)/2
        sq<-(q3-q1)*2
        u<-(qnorm(u,xb,sx)-md)/sq
        return(u)}

qiqe<-function(u,x){ 
	#returns a vector of standardize exponentials (MD and 2IQR)
         #with mean=mean(x) 
        n<-length(x)
        u1<-NULL
        for(i in 1:length(u)){
                if(u[i]>(1/(2*n)) & u[i]<(1-1/(2*n))){
                        u1=c(u1,u[i])
                        }        
        }
                lam<-abs(1/mean(x))
                q1<-qexp(.25,lam)
                q3<-qexp(.75,lam)
                md<-(q3+q1)/2
                sq<-(q3-q1)*2
                xe<-(qexp(u1,lam)-md)/sq
                return(xe)
        
}


qiqp<-function(x){
        data<-x
        x<-sort(qiq(x))
        u<-NULL
        n<-length(x)
        for(i in 1:n) 
                u[i]=(i-.5)/n
        plot(u,x,main="QIQ plot",ylab="data",xlab="P",pch=20,cex=2,
		ylim=c(-1,1))
        abline(-.5,0,col="blue")

        abline(.5,0,col="blue")
        points(c(.25,.75),c(qiq1(x,.25),qiq1(x,.75)),type="l",
		col="red",cex=4)
        u<-NULL
        un<-NULL
        ue<-NULL
        for(i in 1:100){
                u=(i-.5)/100
                if(u>1/(2*n) & u<1-1/(2*n)){
                        un=c(un,u)
                }
        }
        points(un,qiqn(un,data),type="l",col="steelblue4")
        points(c(.02,.05),c(.85,.85),type="l",col="steelblue4")
        text(.07,.85,"normal quantiles",pos=4)
        points(un,qiqe(un,data),type="l",lty=2)
        points(c(.02,.05),c(.95,.95),type="l",lty=2)
        text(.07,.95,"exponential quantiles",pos=4)
        
}


qsig<-function(x){   #Quantile plot of standard deviation (sigma)
        s2<-var(x)
        n<-length(x)
        v<-n-1
        P<-seq(1/(2*n),1-1/(2*n),length=102)
        P<-P[2:100]
        qs<-v*s2/qchisq(1-P,v)
        plot(P,qs,main="Confidence Quantiles of
Sigma",ylab=expression(sigma),type="l",sub="Assuming Normality")
        h<-v*s2/qchisq(.025,v)
        h2<-v*s2/qchisq(.975,v)
        m<-v*s2/qchisq(.5,v)
        md<-(v*s2/qchisq(.25,v)+v*s2/qchisq(.75,v))/2
        h1<-min(qs)-.05*(max(qs)-min(qs))
        points(c(.975,.975,-.1),c(h1,h,h),type="l",col="violetred")
        points(c(.025,.025,-.1),c(h1,h2,h2),type="l",col="violetred")
        points(c(.5,.5,-.1),c(0,m,m),type="l",col="green")
        points(c(1.1,-.1),c(md,md),type="l",col="blue",lty=2)
        points(c(1.1,-.1),c(s2,s2),type="l",col="grey25",lty=2)
        h<-max(qs)-.05*(max(qs)-min(qs))
        legend(.1,.95*max(qs),c("median","MQ","mean"),
		col=c("green","blue","grey25"),lty=c(1,2,2))
}

qmean<-function(x){  #Quantile plot of mean (mu)
        s2<-var(x)
        n<-length(x)
        v<-n-1
        P<-seq(1/(2*n),1-1/(2*n),length=102)
        P<-P[2:100]
        me<-mean(x)
        ub<-me+qt(P,v)*sqrt(s2/n)
        plot(P,ub,main="Confidence Quantiles of
Mean",ylab=expression(mu),type="l",sub="Assuming Normality")
        h<-me+qt(.975,v)*sqrt(s2/n)
        h2<-me+qt(.025,v)*sqrt(s2/n)
        h1<-min(ub)-.05*(max(ub)-min(ub))
        m<-me+qt(.5,v)*sqrt(s2/n)
        ll<-me-(max(x)-min(x))
        md<-(me+qt(.25,v)*sqrt(s2/n)+me+qt(.75,v)*sqrt(s2/n))/2
        points(c(.975,.975,-.1),c(h1,h,h),type="l",col="violetred")
        points(c(.025,.025,-.1),c(h1,h2,h2),type="l",col="violetred")
        points(c(.5,.5,-.1),c(ll,m,m),type="l",col="green")
        points(c(1.1,-.1),c(md,md),type="l",col="blue",lty=2)
        points(c(1.1,-.1),c(me,me),type="l",col="grey25",lty=2)
        h<-max(ub)-.05*(max(ub)-min(ub))
        legend(.1,h,c("median","MQ","mean"),
		col=c("green","blue","grey25"),lty=c(1,2,2))
}


qhist<-function(x,title="quantile histogram"){ #Quantile histogram
        n<-length(x)
        q3<-Qc(x,.75)
        q1<-Qc(x,.25)
        mq<-(q3+q1)/2
        iqr<-(q3-q1)
        endp<-c(mq-2*iqr,mq-1.5*iqr,mq-iqr,mq-.5*iqr,mq,mq+.5*iqr,
		mq+iqr,mq+1.5*iqr,mq+2*iqr)
        bin<-rep(0,10)
        for(i in 2:9){
                for(j in 1:n){
                        if(x[j]<endp[i] && x[j]>endp[i-1])
                                bin[i]=bin[i]+1
                        if(x[j]==endp[i]){
                                bin[i]=bin[i]+.5
                                bin[i+1]=bin[i+1]+.5
                        }
                }
        }

        for(j in 1:n){
                if(x[j]<endp[1]){
                        bin[1]=bin[1]+1
                }
        }
        for(j in 1:n){
                if(x[j]>endp[9]){
                        bin[9]=bin[9]+1
                }
        }
        binp<-bin/n
        plot(c(mq-3*iqr,endp),binp,type="s",
        	xlim=c(min(endp)-.1*(max(endp)-min(endp)),
			max(endp)+.1*(max(endp)-min(endp))),
                xlab="variable",ylab="density",main=title,axes=F)
        points(c(endp[1],endp[1]),c(0,binp[1]),type="l")
        points(c(endp[2],endp[2]),c(0,binp[2]),type="l")
        points(c(endp[3],endp[3]),c(0,binp[3]),type="l")
        points(c(endp[4],endp[4]),c(0,binp[4]),type="l")
        points(c(endp[5],endp[5]),c(0,binp[5]),type="l")
        points(c(endp[6],endp[6]),c(0,binp[6]),type="l")
        points(c(endp[7],endp[7]),c(0,binp[7]),type="l")
        points(c(endp[8],endp[8]),c(0,binp[8]),type="l")
        points(c(endp[9],endp[9]),c(0,binp[9]),type="l")
        points(c(endp[10],endp[10]),c(0,binp[10]),type="l")
        axis(1,rd(c(endp[1],endp[3],endp[5],endp[7],endp[9])))
        axis(2,seq(0,.5,length=11))
}

rqplot<-function(x){  #plots raw quantile plot along with QIQ measurements
        x<-sort(x)
        n<-length(x)
        q3<-Qc(x,.75)
        q1<-Qc(x,.25)
        mq<-(q3+q1)/2
        iqr<-(q3-q1)
        sq<-2*iqr
        ul<-mq+.5*iqr
        ll<-mq-.5*iqr
        l1<-umid(x)
        xl<-NULL
        for(i in 1:n){
                xl[i]=(i-.5)/n
        }
        plot(xl,x,main="raw quantile plot",xlim=c(0,1),
                ylim=c(mq-sq,mq+sq),pch=20,axes=F,ylab="quantile",xlab="P")
        axis(2,c((mq-sq),(mq-iqr),(mq),(mq+iqr),(mq+sq)),
                labels=c(rd(mq-sq),rd(mq-iqr),rd(mq),rd(mq+iqr),rd(mq+sq)))
        axis(1,(0:10)/10,(0:10)/10)
        points(c(.25,.75),c(ul,ul),type="l",col="purple3")
        points(c(.25,.75),c(ll,ll),type="l",col="purple3")
        points(c(.25,.25),c(ul,ll),type="l",col="purple3")
        points(c(.75,.75),c(ll,ul),type="l",col="purple3")
        points(c(.25,.75),c(Qc(x),Qc(x)),type="l",col="maroon4")
        points(c(.25,.75),c(mq,mq),type="l",col="steelblue3",lty=2)
        legend(0,mq+sq,c("median","MQ","Q3 &
Q1"),col=c("maroon4","steelblue3","purple3"),lty=c(1,2,1))
}

procunivariate<-function(x,name="data",type=1){
        if(type==1){
                op<-par(mfrow=c(1,3),mar=c(5,4,4,1),oma=c(0,0,3,0),
			cex.main=1.5,font.main=1)
                rqplot(x)
                qiqp(x)
                qhist(x)
                title(name,outer=T,cex.main=3,font.main=6)
                ck <- readline('strike "Enter" to continue')
                op<-par(mfrow=c(1,2),cex.main=1,font.main=1,mar=c(6,4,4,1))
                qmean(x)
                qsig(x)
                title(name,outer=T,cex.main=2,font.main=6)
                }
        if(type==2){
                op<-par(mfrow=c(2,3),mar=c(5,4,4,1),
			cex.main=1.5,font.main=1)
                qufplot(x)
                qiqp(x)
                qsig(x)
                qmean(x)
                qhist(x)
                rqplot(x)}
ck <- readline('strike a "Enter" to continue')
n<-length(x)
q3<-quantile(x,P_(.75,n))
q1<-quantile(x,P_(.25,n))
MQ<-(q3+q1)/2
IQR<-q3-q1
M<-summary(x)
M1<-M
M1[4]=M[5]
M1[5]=M[6]
M1[10]=M[4]
M=M1
M[6]=qiq1(x,.5/length(x))
M[7]=qiq1(x,.5)
M[8]=qiq1(x,(length(x)-.5)/length(x))
M[9]=MQ
M[10]=IQR*2
M[11]=mean(x)
M[12]=sqrt(var(x))
lis<-c("Min =","1st Qt =","Med =","3rd Qt =","Max =","QI(Min) =",
		"QI(.5) =","QI(Max)=","MQ =","2IQR =","mean","std")
names(M)=lis
op<-par(mfrow=c(1,1),mar=c(1,1,4,1),cex.main=1.5,font.main=1)
plot(1:100,1:100,type="n",axes=F,main="EQA summary",xlab="",ylab="")
text(rep(25,5),seq(90,10,by=-20),lis[1:5],pos=2)
text(rep(25,5),seq(90,10,by=-20),round(M[1:5],5),pos=4)
text(rep(75,5),seq(90,10,by=-20),lis[6:10],pos=2)
text(rep(75,5),seq(90,10,by=-20),round(M[6:10],5),pos=4)

ck <- readline('strike "Enter" to continue')

op<-par(mfrow=c(1,1),mar=c(1,1,4,1),cex=1.5,oma=c(0,0,0,0))
plot(1:100,1:100,type="n",axes=F,main="Quantile Points (Normal)",xlab="",
	ylab="")
pm<-c("p",.5,.25,.1,.05,.025,.01,.005)
pl<-c(.5,.25,.1,.05,.025,.01,.005)
q1<-NULL
q2<-NULL
sq1<-NULL
sq2<-NULL
for(i in 1:7){
        if(pl[i]>1/(2*n) & pl[i]<1-1/(2*n)){  
		#checking for quantile estimability
                q1[i]=rd(qnorm(pl[i]),4)
                q2[i]=rd(qnorm(1-pl[i]),4)
                sq1[i]=rd(Qc(x,pl[i]),4)
                sq2[i]=rd(Qc(x,1-pl[i]),4)
        }
        else{
                q1[i]="------"
                q2[i]="------"
                sq1[i]="------"
                sq2[i]="------"
        }
}

text(rep(15,5),seq(90,10,length=8),pm,pos=2)
op<-par(cex=1)
text(rep(30,5),seq(90,10,length=8),c("Q(p;N(0,1))",q1),pos=4)
text(rep(45,5),seq(90,10,length=8),c("Q(p;data)",sq1),pos=4)
text(rep(70,5),seq(90,10,length=8),c("Q(1-p;N(0,1))",q2),pos=4)
text(rep(85,5),seq(90,10,length=8),c("Q(1-p; data)",sq2),pos=4)
 points(c(23,23),c(0,100),type="l")
abline(85,0)


ck <- readline('strike "Enter" to continue')
op<-par(cex=1.5)
pm<-c("p",.05,.25,.5,.75,.95)
pl<-c(.05,.25,.5,.75,.95)
plot(1:100,1:100,type="n",axes=F,main="Quantile Points (Exponential)",
	xlab="",ylab="")
q1<-NULL
sq1<-NULL
for(i in 1:5){
        if(pl[i]>1/(2*n) & pl[i]<1-1/(2*n)){  
	#checking for quantile estimability
        q1[i]=rd(qexp(pl[i]),4)
        sq1[i]=rd(Qc(x,pl[i]),4)
        }
        else{
        q1[i]="------"
        sq1[i]="------"
        }
}
text(rep(20,5),seq(90,10,length=6),pm,pos=2)
text(rep(40,5),seq(90,10,length=6),c("Q(p;Exp(1))",q1),pos=4)
text(rep(70,5),seq(90,10,length=6),c("Q(p;data)",sq1),pos=4)
points(c(30,30),c(0,100),type="l")
abline(85,0)
}

procunivariate(rnorm(30),"30 Random Normals")