###############################################################################
####This program computes all the results in Section 2 about distribution  ####
####free tolerance/reference rectangles by using Tukey's equivalence blocks####
####The sample size nn must be large enough in order for the tolerance     ####
####rectangle to be constructed.                                           ####
####nn: the sample size;                                                   ####
####rp: the content;                                                       ####
####rgamma: the confidence level;                                           ####
####xx: the sample available.                                              ####
###############################################################################


##Plot Figure 1##
#################
xx1 <- c(5, 3,  1,  2, 4, 2.5, 3.5)
xx2 <- c(3, 5, 3.5, 1, 2, 4,    1.5)
plot(xx1,xx2,pch=4, col="red", xlim=c(0, 6),
     ylim=c(0, 6), xlab = expression(x[1]),ylab = expression(x[2]))
#plot(xx1,xx2,pch=4, col="red", xlim=c(0, 6),
#     ylim=c(0, 6), xlab = tex('$x_1$'),ylab = tex('$x_2$'))
lines(c(5,5), c(0,6), col="blue")
text(5.5, 3, expression(B1))
lines(c(0,5), c(5,5), col="blue")
text(2.5, 5.5, expression(B2))
lines(c(1,1), c(0,5), col="blue")
text(0.3, 3, expression(B3))
lines(c(1,5), c(1,1), col="blue")
text(3, 0.5, expression(B4))
lines(c(4,4), c(1,5), col="blue")
text(4.5, 3, expression(B5))
lines(c(1,4), c(4,4), col="blue")
text(2.5, 4.5, expression(B6))
lines(c(3.5,3.5), c(1,4), col="blue")
text(2.25, 2.5, expression(B7))
text(3.75, 2.5, expression(B8))


##Define the function that finds the smallest sample size nn       ##
##so that a (rp, rgamma) tolerance rectangle can be constructed via##
##the union of several equivalence blocks                          ##
##nk2: minimum sample size for the 2-sided tolerance rectangle     ##
##nk1: minimum sample size for the 1-sided tolerance rectangle     ##
#####################################################################
minSampSize <- function(rp, rgamma){
  n0 = ceiling( log(1-rgamma)/log(rp) )

  nk1=max(n0,2)
  while ( pbeta(rp, nk1-1, 2, lower.tail=FALSE) < rgamma ){
    nk1 <- nk1 + 1
  }  #the minimum sample size required for 1-sided
  
  nk2=max(n0,4)
  while ( pbeta(rp, nk2-3, 4, lower.tail=FALSE) < rgamma ){
    nk2 <- nk2 + 1
  }  #the minimum sample size required for 2-sided
  
  return(list("min sample size for 1-s",nk1,"min sample size for 2-s",nk2))
}
rp=0.95 #rp=0.90  rp=0.99
rgamma = 0.95  #rgamma=0.90  rgamma=0.99 
minSampSize(rp, rgamma)
pbeta(rp, minSampSize(rp, rgamma)[[2]][1]-1, 2, lower.tail=FALSE) #this should be >= rgamma
pbeta(rp, minSampSize(rp, rgamma)[[2]][1]-2, 2, lower.tail=FALSE) #this should be < rgamma
pbeta(rp, minSampSize(rp, rgamma)[[4]][1]-3, 4, lower.tail=FALSE) #this should be >= rgamma
pbeta(rp, minSampSize(rp, rgamma)[[4]][1]-4, 4, lower.tail=FALSE) #this should be < rgamma


##Define the function that finds k0 so that the union of k0##
##equivalence blocks is a (rp, rgamma) tolerance region    ##
#############################################################
numEquiBloc<-function(rp,rgamma,nn){
  if(nn<min(minSampSize(rp,rgamma)[[4]][1], minSampSize(rp,rgamma)[[2]][1])){
    print( "The sample size is too small even for 1-sided rectangle")
    return("The sample size is too small even for 1-sided rectangle")
  }else{
    nk=nn
    while ( pbeta(rp, nk, nn+1-nk, lower.tail=FALSE) > rgamma ){
      nk <- nk - 1
    }
    k = nk + 1   #The  k  required
    return(k)
  } 
}
nn = 200
rp = 0.95
rgamma=0.95
minSampSize(rp, rgamma)
nk=numEquiBloc(rp,rgamma,nn)
pbeta(rp, nk, nn+1-nk, lower.tail=FALSE)   #this is the true confLevel, sould be > rgamma
pbeta(rp, nk-1, nn+2-nk, lower.tail=FALSE) #this sould be < rgamma

#For the simple example with n=7 given at the beginning
nn=7
rp=0.45 
rgamma=0.60
minSampSize(rp, rgamma)
nk=numEquiBloc(rp,rgamma,nn)
pbeta(rp, nk, nn+1-nk, lower.tail=FALSE)   #this is the true confLevel, sould be > rgamma
pbeta(rp, nk-1, nn+2-nk, lower.tail=FALSE) #this sould be < rgamma


#For the discussion in the last paragraph of Section 2 with n=100
nn=1000
rp=0.90 
rgamma=0.90
minSampSize(rp, rgamma)
nk=numEquiBloc(rp,rgamma,nn)
pbeta(rp, nk, nn+1-nk, lower.tail=FALSE)   #this is the true confLevel, sould be > rgamma
pbeta(rp, nk-1, nn+2-nk, lower.tail=FALSE) #this sould be < rgamma


#######################################################################
##Define the main function that finds the (rp,rgamma) tolerance region##
##by deleting nn+1-nk equivalence blocks, in the order of R,T,L,B.   ##
##rp: the content;                                                   ##
##rgamma: the confidence level;                                       ##
##xx: the sample available in the form of a nn*2 matrix.             ##
#######################################################################
RectRTLB2 <- function(xx,rp,rgamma){
  nn <- dim(xx)[1]
  nk <- numEquiBloc(rp,rgamma,nn)
  nkD <- nn+1-nk #number of equiv blocks to delete
  nkDCycl <- floor(nkD/4) #number of cycles for striping 4 blocks R,T,L,B
  nkDRem <- nkD - 4*nkDCycl #number of blocks (<4) to strip further
  
  for (i in 1:nkDCycl) { #Strip away nkDCycl cycles of  4 blocks R,T,L,B

    maxRind <- which(xx[,1]==max(xx[,1])) #Indices of the obs having the max x_1 value
    maxRval <- xx[maxRind[1],1] #Record the max x_1 value
    xx <- xx[-maxRind[1],]  #Delete this observation from the sample
    
    maxTind <- which(xx[,2]==max(xx[,2])) #Indices of the obs having the max x_2 value
    maxTval <- xx[maxTind[1],2] #Record the max x_2 value
    xx <- xx[-maxTind[1],]  #Delete this observation from the sample
    
    minLind <- which(xx[,1]==min(xx[,1])) #Indices of the obs having the min x_1 value
    minLval <- xx[minLind[1],1] #Record the min x_1 value
    xx <- xx[-minLind[1],]  #Delete this observation from the sample
    
    minBind <- which(xx[,2]==min(xx[,2])) #Indices of the obs having the min x_2 value
    minBval <- xx[minBind[1],2] #Record the min x_2 value
    xx <- xx[-minBind[1],]  #Delete this observation from the sample
  }
  
  if(nkDRem == 0){ #Strp the remaining (<4) block(s)
    "Do nothing"
  }else if(nkDRem == 1){
    maxRind <- which(xx[,1]==max(xx[,1])) #Indices of the obs having the max x_1 value
    maxRval <- xx[maxRind[1],1] #Record the max x_1 value
    xx <- xx[-maxRind[1],]  #Delete this observation from the sample
  }else if(nkDRem == 2){
    maxRind <- which(xx[,1]==max(xx[,1])) #Indices of the obs having the max x_1 value
    maxRval <- xx[maxRind[1],1] #Record the max x_1 value
    xx <- xx[-maxRind[1],]  #Delete this observation from the sample
    
    maxTind <- which(xx[,2]==max(xx[,2])) #Indices of the obs having the max x_2 value
    maxTval <- xx[maxTind[1],2] #Record the max x_2 value
    xx <- xx[-maxTind[1],]  #Delete this observation from the sample
  }else if(nkDRem == 3){
    maxRind <- which(xx[,1]==max(xx[,1])) #Indices of the obs having the max x_1 value
    maxRval <- xx[maxRind[1],1] #Record the max x_1 value
    xx <- xx[-maxRind[1],]  #Delete this observation from the sample
    
    maxTind <- which(xx[,2]==max(xx[,2])) #Indices of the obs having the max x_2 value
    maxTval <- xx[maxTind[1],2] #Record the max x_2 value
    xx <- xx[-maxTind[1],]  #Delete this observation from the sample
    
    minLind <- which(xx[,1]==min(xx[,1])) #Indices of the obs having the min x_1 value
    minLval <- xx[minLind[1],1] #Record the min x_1 value
    xx <- xx[-minLind[1],]  #Delete this observation from the sample
  }    
  return(list( c(maxRval,maxTval,minLval,minBval), c(nk, nkD) ))
}

##Apply to one simulated dataset as an example##
set.seed(5)
nn=250
xx <- matrix(rnorm(nn*2),ncol=2)
rp=0.95
rgamma=0.95
bounds<-RectRTLB2(xx,rp,rgamma)[[1]]  #The (R, T, L, B) values of the tolerance region

plot(xx[,1],xx[,2],pch=4, col="red", xlim=c(min(xx[,1])-0.5, max(xx[,1])+0.5),
      ylim=c(min(xx[,2])-0.5, max(xx[,2])+0.5), xlab = expression(x[1]),ylab = expression(x[2]))
xxx <- c(bounds[3],bounds[1],bounds[1],bounds[3])
yyy <- c(bounds[4],bounds[4],bounds[2],bounds[2])
polygon(xxx, yyy, col = NA, border = "blue")
title("Tolerance rectangle")


##Plot Figure 2##
#################
xx1 <- c(5, 3,  1,  2, 4, 2.5, 3.5)
xx2 <- c(3, 5, 3.5, 1, 2, 4,    1.5)
plot(xx1,xx2,pch=4, col="red", xlim=c(0, 6),
     ylim=c(0, 6), xlab = expression(x[1]),ylab = expression(x[2]))
#plot(xx1,xx2,pch=4, col="red", xlim=c(0, 6),
#     ylim=c(0, 6), xlab = tex('$x_1$'),ylab = tex('$x_2$'))
lines(c(5,5), c(0,6), col="blue")
text(5.5, 3, expression(B1))
lines(c(0,5), c(5,5), col="blue")
text(2.5, 5.5, expression(B2))
lines(c(4,4), c(0,5), col="blue")
text(4.5, 2.5, expression(B3))
lines(c(0,4), c(4,4), col="blue")
text(2, 4.5, expression(B4))
lines(c(3.5,3.5), c(0,4), col="blue")
text(3.75, 2, expression(B5))
lines(c(0,3.5), c(3.5,3.5), col="blue")
text(1.75, 3.75, expression(B6))
lines(c(2,2), c(0,3.5), col="blue")
text(2.75, 1.75, expression(B7))
text(1, 1.75, expression(B8))


rp = 0.95
rgamma=0.95
nn <- minSampSize(rp, rgamma)[[2]][1]
nk=numEquiBloc(rp,rgamma,nn)
pbeta(rp, nk, nn+1-nk, lower.tail=FALSE)   #this is the true confLevel, sould be > rgamma
pbeta(rp, nk-1, nn+2-nk, lower.tail=FALSE) #this sould be < rgamma

rp = 0.45
rgamma=0.60
minSampSize(rp, rgamma)[[2]][1]
nn = 7
nk=numEquiBloc(rp,rgamma,nn)
pbeta(rp, nk, nn+1-nk, lower.tail=FALSE)   #this is the true confLevel, sould be > rgamma
pbeta(rp, nk-1, nn+2-nk, lower.tail=FALSE) #this sould be < rgamma


##########################################################################
##Define the main function that finds the 1-sided (rp,rgamma) tolerance  ##
##rectangle by deleting nn+1-nk equivalence blocks, in the order of R,T.##
##rp: the content;                                                      ##
##rgamma: the confidence level;                                          ##
##xx: the sample available in the form of a nn*2 matrix.                ##
##########################################################################
RectRTLB1 <- function(xx,rp,rgamma){
  nn <- dim(xx)[1]
  nk <- numEquiBloc(rp,rgamma,nn)
  nkD <- nn+1-nk #number of equiv blocks to delete
  nkDCycl <- floor(nkD/2) #number of cycles for striping 2 blocks R,T
  nkDRem <- nkD - 2*nkDCycl #number of blocks (<2) to strip further
  
  for (i in 1:nkDCycl) { #Strip away nkDCycl cycles of  2 blocks R,T
    
    maxRind <- which(xx[,1]==max(xx[,1])) #Indices of the obs having the max x_1 value
    maxRval <- xx[maxRind[1],1] #Record the max x_1 value
    xx <- xx[-maxRind[1],]  #Delete this observation from the sample
    
    maxTind <- which(xx[,2]==max(xx[,2])) #Indices of the obs having the max x_2 value
    maxTval <- xx[maxTind[1],2] #Record the max x_2 value
    xx <- xx[-maxTind[1],]  #Delete this observation from the sample
  }
  
  if(nkDRem == 0){ #Strp the remaining (<2) block(s)
    "Do nothing"
  }else if(nkDRem == 1){
    maxRind <- which(xx[,1]==max(xx[,1])) #Indices of the obs having the max x_1 value
    maxRval <- xx[maxRind[1],1] #Record the max x_1 value
    xx <- xx[-maxRind[1],]  #Delete this observation from the sample
  }    
  return(list( c(maxRval,maxTval), c(nk, nkD) ))
}

##Apply to one simulated dataset as an example##
set.seed(5)
nn=250
xx <- matrix(rnorm(nn*2),ncol=2)
rp=0.95
rgamma=0.95
bounds<-RectRTLB1(xx,rp,rgamma)[[1]]  #The (R, T, L, B) values of the tolerance region

plot(xx[,1],xx[,2],pch=4, col="red", xlim=c(min(xx[,1])-0.5, max(xx[,1])+0.5),
     ylim=c(min(xx[,2])-0.5, max(xx[,2])+0.5), xlab = expression(x_1),ylab = expression(x_2))
lines(c(-4,bounds[1]),c(bounds[2],bounds[2]), col="blue")
lines(c(bounds[1],bounds[1]),c(-3.4,bounds[2]), col="blue")
title("One-sided tolerance rectangle")


##Define the function that finds the smallest sample size nn     ##
##so that a (rp, rgamma) tolerance boxes can be constructed via  ##
##the union of several equivalence blocks                        ##
##nk2: minimum sample size for the 2-sided tolerance box         ##
##nk1: minimum sample size for the 1-sided tolerance box         ##
###################################################################
minSampSize3 <- function(rp, rgamma){   #for p=3 dimensions only
  n0 = ceiling( log(1-rgamma)/log(rp) )
  
  nk1=n0
  while ( pbeta(rp, nk1-2, 3, lower.tail=FALSE) < rgamma ){
    nk1 <- nk1 + 1
  }  #the minimum sample size required for 1-sided
  
  nk2=n0
  while ( pbeta(rp, nk2-5, 6, lower.tail=FALSE) < rgamma ){
    nk2 <- nk2 + 1
  }  #the minimum sample size required for 2-sided
  
  return(list("min sample size for 1-s",nk1,"min sample size for 2-s",nk2))
}
rp=0.95  #rp=0.90  rp=0.99
rgamma = 0.95  #rgamma = 0.90  rgamma = 0.99
minSampSize3(rp, rgamma)
pbeta(rp, minSampSize3(rp, rgamma)[[2]][1]-2, 3, lower.tail=FALSE) #this should be >= rgamma
pbeta(rp, minSampSize3(rp, rgamma)[[2]][1]-3, 3, lower.tail=FALSE) #this should be < rgamma
pbeta(rp, minSampSize3(rp, rgamma)[[4]][1]-5, 6, lower.tail=FALSE) #this should be >= rgamma
pbeta(rp, minSampSize3(rp, rgamma)[[4]][1]-6, 6, lower.tail=FALSE) #this should be < rgamma