################################################################################# ####This program computes the critical constant rc for the bivariate-normal #### ####p-content 1-alpha confidence tolerance ellipse, by using the new method, #### ####the KM method available in the package tolerance and the MM method, #### ####for given rp-content, rconf-confidence level and nn-sample size. #### #### ################################################################################# ####Clear the workspace first#### ls() rm(list=ls()) pwd() ################################################################################ ##The prob in exps (6) minus the content, i.e. the function to solve root from## ################################################################################ funProb2 <- function(cc, vW, vlam, rp){ # cc = c/(n-1) ##Define the two integration limits## ##################################### xL <- vW[1] - sqrt( cc/vlam[1] ) #the lower limit of v1 xU <- vW[1] + sqrt( cc/vlam[1] ) #the upper limit of v1 ##Define the integrand## ######################## funIntegrand <- function(v1){ delt <- sqrt( (cc-vlam[1]*(v1-vW[1])^2)/vlam[2] ) #works y <- dnorm(v1)* ( pnorm(vW[2]+delt) - pnorm(vW[2]-delt) ) return(y) } y <- suppressWarnings( integrate(funIntegrand, lower=xL, upper=xU)[[1]] ) - rp ##[[1]] is essential to get the value of the integration return(y) } #funProb2(cc, vW, vlam, rp) ########################################################################### ###Define the function that finds the critical constant rc by using ### ### the new method proposed. ### ###rp: the content of the tolerance ellipsoid ### ###rconf: the confidence level of the tolerance ellipsoid ### ###nn: the sample size of the sample available ### ###nsimu: number of simulations ### ###tol: error tolerance in the bisection method ### ###Nseed: the seed number used in set.seed ### ########################################################################### Frc <- function(rp,rconf,nn,nsimu, tol, Nseed){ set.seed(Nseed) rnm1 <- 1/(nn-1) rsqrtn <- 1/sqrt(nn) rc = rep(0,nsimu) #Compute all the c's for each of nsimu simulation replications for (i in 1:nsimu) { #mLamd <- rNonsingularWishart(1, nn-1, diag(c(1,1)) )[,,1] #Generate \Lambda c1 <- sqrt(rchisq(1,df=nn-1)) c2 <- sqrt(rchisq(1,df=nn-2)) nij <- rnorm(1) mG <- matrix(c(c1,0,nij,c2),nrow=2,ncol=2) mLamd <- t(mG)%*%mG mM <- solve(mLamd) #The matrix M which is the inverse of \Lambda #mM%*%mLamd #Should be equal to the identity matrix vZ <- matrix(rnorm(2), ncol = 1) #The column vector z ##Find the c for this simulation replication using the bisection method## #Orthogonal diagonization of matrix mM resu <- eigen(mM, symmetric=TRUE) #Find the decomposition of matrix mM mQ <- t(resu$vectors) #the matrix Q vlam <- resu$values #the eigen values #t(mQ)%*%diag(vlam)%*%mQ #this should be equal to mM #t(mQ)%*%mQ #this should be equal to the identity matrix, i.e. mQ is orthogonal #mQ%*%mM%*%t(mQ) #this should be equal to diag(vlam) vW <- rsqrtn* mQ%*%vZ #the input for funProb2 #cc <- c * rnm1 #First find the interval (clow, cupp) that contains cc <- c * rnm1 Qchi <- qchisq( rp, df=2, ncp = t(vZ)%*%vZ /nn ) #clow <- (nn-1)*min(vlam)*Qchi #cupp <- (nn-1)*max(valm)*Qchi clow <- min(vlam)*Qchi cupp <- max(vlam)*Qchi flow <- funProb2(clow, vW, vlam,rp) #Make sure this < 0 fupp <- funProb2(cupp, vW, vlam,rp) #Make sure this > 0 cmid <- (clow+cupp)/2 while ( (cupp-clow)>tol ){ cmid <- (clow+cupp)/2 c <- cmid fmid <- funProb2(c, vW, vlam,rp) if ( fmid == 0 ){ rc[i] = cmid break } else if (flow*fmid <0){ cupp <- cmid fupp <- fmid } else { clow <- cmid flow <- fmid } } rc[i] = cmid #cmid } #cc = quantile(rc, rconf) cc <- (nn-1) * sort(rc)[nsimu*rconf] # This is the critical constant sought return(cc) } nsimu <- 1000000 #Number of simulations tol <- 0.0001 #absolute error of the solution c Nseed <- 5 rp <- 0.90 #content covered for example 1 of Section 4 rp <- 0.95 #content covered for example 1 of MM (2015) in Section 4 rconf <- 0.95 #confidence level nn <- 30 #sample size pc<-proc.time() # Keep track of computation time rc_NEW <- Frc(rp,rconf,nn,nsimu, tol, Nseed) ctime <- proc.time()-pc c(rc_NEW, ctime[1]) ####################################################################### ###Compute the critical constant rc by using the R package tolerance### ### -- the KM method, which is recommended in tolerance ### ####################################################################### x <- matrix( rep(1,nn*2), ncol = 2) library(tolerance) set.seed(5) pc<-proc.time() out1 <- mvtol.region(x = x, alpha = 1-rconf, P = rp, B = nsimu, method = "KM") rc_KM <- out1[[1]] ctime <- proc.time()-pc time_KM <- ctime[1] c(rc_KM, time_KM) ############################################################### ###Compute the critical constant rc by using the MM method ### ### ############################################################### ##Define a function that is used in the main function of the MM method## ######################################################################## funADF <- function(c){ #the values of (a, delta, f) for given c=(c1,c2,c3,c4) s1 <- c[3]/(c[2])^(3/2) s2 <- c[4]/(c[2])^2 if(s1^2 > s2){ ra <- 1/(s1-sqrt(s1^2-s2)) rdelta <- s1*ra^3 - ra^2 rf <- ra^2 - 2*rdelta }else{ ra <- 1/s1 rdelta <- 0 rf <- ra^2 } y <- c(ra,rdelta,rf) return(y) } ###Define the function that computes the critical constant rc by MM method### ############################################################################# rcMM <- function(rp,rconf,nn,nsimu,Nseed){ set.seed(5) #Generate all the simulation data once SimDat <- matrix(0, nrow=5,ncol = nsimu) #One column is used in one simulation repetition SimDat[1,] <- rchisq(nsimu,df=nn-1) #1st three elements in each column is for generating M = Lambda^{-1} SimDat[2,] <- rchisq(nsimu,df=nn-2) SimDat[3,] <- rnorm(nsimu) SimDat[4,] <- rnorm(nsimu) #The last two elements in each column is for generating \bz SimDat[5,] <- rnorm(nsimu) rc = rep(0,nsimu) #Compute all the c's for each simulation replication for (i in 1:nsimu) { m11 = ( SimDat[2,i] + (SimDat[3,i])^2 )/(SimDat[1,i] * SimDat[2,i]) #find the matrix M in (8) m12 = - sqrt(SimDat[1,i]) * SimDat[3,i] / (SimDat[1,i] * SimDat[2,i]) m22 = SimDat[1,i] / (SimDat[1,i] * SimDat[2,i]) mM <- matrix(c(m11,m12,m12,m22), nrow=2) resu <- eigen(mM, symmetric=TRUE) #Find the decomposition of mM V <- resu$vectors #the matrix L lam <- resu$values #the eigen values #V%*%diag(lam)%*%t(V) #this should be equal to mM #V%*%t(V) #this should be equal to the identity matrix #t(V)%*%mM%*%V #this should be equal to diag(lam) z1 <- SimDat[4,i] # the z1 and z2 in (9) and (10) z2 <- SimDat[5,i] vdelta <- (t(V)%*%matrix(c(z1,z2), nrow=2))^2 / nn # the vector of non-centrality parameters c(delta_i) ccc = rep(0,4) # The c_j values in this simulation replication for(j in 1:4){ ccc[j] <- sum(lam^j) + j*sum( (lam^j)*vdelta ) } ADF <- funADF(ccc) # the values of (a, delta, f) rc[i] = (nn-1)*( sqrt(ccc[2]/((ADF[1])^2)) * ( qchisq(rp, df=ADF[3], ncp = ADF[2]) - ADF[3]- ADF[2] ) + ccc[1] ) #Using the formula } #cc = quantile(rc, rconf) cc <- sort(rc)[nsimu*rconf] # This is the critical constant sought return(cc) } rp <- 0.90 #content covered rconf <- 0.95 #confidence level nn <- 30 #sample size nsimu <- 100000 #Number of simulations Nseed <- 5 pc<-proc.time() # Keep track of computation time rc_MM <- rcMM(rp,rconf,nn,nsimu,Nseed) ctime <- proc.time()-pc time_MM <- ctime[1] c(rc_MM, time_MM) print(c(rc_NEW, rc_KM, rc_MM)) #write(c(rc_NEW, rc_KM, rc_MM), file="Outcome.txt", ncolumns=3, append=TRUE) #cat("at iteration",i,"the c is",out1,"\n")