Abstract tubes associated with perturbed polyhedra with applications to multidimensional normal probability computations

1 Abstract tub es asso ciated with p erturb ed p olyhedra with applications to m ultidime nsional normal probability computations S. KURIKI The Institute of Statistic al Mathematics 10-3 Midoricho, T achikawa, T okyo 190-8562, Jap an E-mail: kuriki@ism.ac.jp T. MIW A National Institute for A gr o-Envir onmental Scie nc es 3-1-3 Kannondai, Tsukub a 305-8604, Jap an E-mail: miwa@niaes.affr c.go.jp A. J. HA YTER Daniels Col le ge of Business, University of De nver 2101 S. University Blvd., Denver, CO 8020 8-8921, USA E-mail: An thony.Hayter@du.e du Let K b e a closed conv ex p olyhedron defined b y a finite num ber of li near in- equalities. In this pap er we refine the theory of abstract tub es (Naiman and Wynn, 1997) asso ciated with K when K is perturb ed. In particular, we fo- cus on the p erturbation that is lexicographic and in an outer di rection. An algorithm for constructing the abstract tube by means of li near programming and its implemen tation are discussed. Using t he abstract tu be for perturb ed K com bined with the recursive inte gration tec hnique proposed by Mi wa, Hayter and Kuriki (2003), we show that the multidimensional normal probability for a polyhedral region K can be computed efficiently . In addition, abstract tubes and the distribution functions of studen tized range statistics are exhibited as n umerical examples. Keywor ds : Abstract T ube; Inclusion-Exclusion Iden tity; Lexicographic Method; Linear Programming; Multiple Comparisons; Perturbation; Studen- tized Range. 1. In tro ductio n Let A = ( a 1 , . . . , a m ) be an n × m matrix such that a i 6 = 0 for a ll i , and let b = ( b 1 , . . . , b m ) ⊤ ∈ R m be an m × 1 constant vector. Define a closed 2 conv ex po ly hedron by K = { x ∈ R n | A ⊤ x ≤ b } . Throughout the pa per, K is assumed to b e no nempt y . Suppose tha t x ∈ R n is a random vector dis tributed as an n -dimensional standar d normal distri- bution N n (0 , I n ). The primar y mo tiv ation and hence one o f the purposes of this study is to ev a luate the n -dimensio nal normal probability P ( K ) = P n ( K ) := Pr( x ∈ K ) . (1) When m = n and A is no n-singular, K is a co ne referr ed to as the simple cone formed as the intersection of n half spaces loca ted in the general po sitions (Ba rvinok, 2 002 1 ). In this pap er, we do not r e q uire that the ap ex of cone lies at the o rigin. In the ca se where m < n a nd rank( A ) = m , K is decomp osed as K = K 1 ⊕ R n − m with K 1 an m -dimensional simple cone, where “ ⊕ ” means the or thogonal direct sum. In this case the probability (1) is obtained as P n ( K ) = P m ( K 1 ). F or such a simple c o ne K , Miw a, Hayter and Kuriki (20 0 3) 2 prop osed a recur siv e integration algor ithm to ev aluate the probability P ( K ) by generaliz ing the idea of Abrahamson (1964). 3 This algorithm is practically useful for dimensions n up to 20, and is av aila ble in the R librar y mvtn orm . 4,5 In this paper, using this integration a lgorithm as a building block, w e demonstrate that the pr o babilit y P ( K ) for any p olyhedro n K can b e com- puted b y means o f the abstra ct tube ideas pr opo sed by Naiman and Wynn (1997). 6 This metho d provides a sophistica ted v er sion o f the inclusion- exclusion identit y , and divide s the complimentary set K c int o sig ned co nes. By calculating the norma l pr obability for each co ne and summing them up, we obtain P ( K ) = 1 − P ( K c ). How ever, as we shall see later , the resulting cones ar e no t necessarily simple, and w e need to in tro duce a p erturbation int o the linear inequa lit y s ystem. W e fo cus on a particula r type of p ertur- bation, and refine the theory of abstract tubes for a p erturb ed inequality system in this setting. Some statements ar e s imply taken from Naiman and Wynn (1 997), 6 but some o f them are nov el. F or example, we show that under our p erturbation the p erturb ed system alwa ys defines an abstract tube instead of a “ w eak abstra ct tub e”. Moreov er, we show that the abstract tub e for the pertur bed system ca n be constructed by mea ns of a standard linear programming technique. W e prop ose that such an a lgorithm should b e used, and we make some re- marks that should b e inco rpo rated into its implemen tation. The prototype R pr ogram is av ailable from the authors. 3 The construction of the pap er is a s follows. In Section 2, the theo ry of the abstract tube for a p erturb ed system is inv estig ated. In Section 3, the algo rithm and its implementation for co nstructing the abstract tube is discussed. In Section 4, we exhibit numerical examples. F or T ukey’s studen- tized r a nge statistics, the abstra ct tub es are co ns tructed in v arious settings , and their distribution functions are computed. 2. Abstract tub es for a p erturbe d inequalit y system In this section we in vestigate the abstract tub e asso ciated with a clo sed conv ex p olyhedron. Because of the reason stated la ter, w e fo cus on the case where the inequality system defining the p olyhedron is pe rturbe d. W e prop ose a lex ic ographic p erturbation in an outer directio n (see (4)) and refine the theoretica l results of Naiman and Wynn (1 997) 6 in this setting. Define m half spaces H i = { x ∈ R n | a ⊤ i x ≤ b i } , i = 1 , . . . , m. The complement a nd the bounda ry of H i are denoted by H c i = { x ∈ R n | a ⊤ i x > b i } and ∂ H i = { x ∈ R n | a ⊤ i x = b i } , resp ectively . Since P ( K ) = P ( ∩ m i =1 H i ) = 1 − P ( ∪ m i =1 H c i ), we ca n ev aluate P ( ∪ m i =1 H c i ) ins tea d of P ( ∩ m i =1 H i ). The family o f the complement s of the half spaces is denoted by H = { H c i | i = 1 , . . . , m } . Let F i = ∂ H i ∩ K , i = 1 , . . . , m. If F i 6 = ∅ , F i is a face of K . Define a family o f subsets of the indices { 1 , . . . , m } by F = { J ⊆ { 1 , . . . , m } | J 6 = ∅ , ∩ i ∈ J F i 6 = ∅} . Then, F forms an (abstract) simplicial complex . {∩ i ∈ J F i | J ∈ F } is the family o f a ll faces (vertex, edge, ..., fac et) of K . Note that ∩ i ∈ J 1 F i = ∩ i ∈ J 2 F i can ha pp en even though J 1 6 = J 2 . Let 1 S ( x ) = ( 1 , if x ∈ S , 0 , otherwise , be the indicator function o f the set S . The following is Theor em 2 of Naiman and Wynn (1997 ). 6 4 Prop osition 2.1. The p air ( H , F ) is an abstr act tub e in the sense t hat 1 ∪ m i =1 H c i ( x ) = X J ∈F ( − 1) | J |− 1 1 ∩ i ∈ J H c i ( x ) for al l x ∈ R n , (2) wher e | J | is the c ar dinality of the set J . Assuming that x ∈ R n is a standard Gaussian vector distributed as N n (0 , I n ), a nd taking the exp ectation of (2), we obtain the formula 1 − P ( K ) = X J ∈F ( − 1) | J |− 1 P ( ∩ i ∈ J H c i ) . In the right side, if ∩ i ∈ J H c i forms an n -dimensional simple cone or the direct sum of a simple cone a nd a linea r subspace, we can use the algorithm of Miwa et al. (2003) 2 for ev aluating the probability P ( ∩ i ∈ J H c i ). Ho wev er , this is not alwa ys the case, a nd for it to b e the ca se the condition of “genera l po sition” is required. B e fo re defining this notion, we e x plain tw o examples of abstract tub es that do not meet this co ndition. Example 2.1. Consider a p olyhedra l cone called a “pyramid” K ⊂ R 3 consisting of ( x 1 , x 2 , x 3 ) s atisfying − x 1 − x 2 + x 3 ≤ 1 [1] − x 1 + x 2 + x 3 ≤ 1 [2] + x 1 + x 2 + x 3 ≤ 1 [3] + x 1 − x 2 + x 3 ≤ 1 [4] The boundar ies of these half spac e s are not in the general p osition in the sense that 4 hyperplane s share a p oint (0 , 0 , 1) in R 3 . The simplicial complex is shown to b e F = { 1 , 2 , 3 , 4 , 1 2 , 13 , 14 , 2 3 , 24 , 3 4 , 123 , 124 , 134 , 23 4 , 1234 } . (3) In the expressio n above “ 134” means { 1 , 3 , 4 } for example, and this c on- ven tion is adopted throughout the pa p er. Applying this to the P rop osition 2.1, we can decomp ose 1 − P ( K ) into |F | = 1 5 terms . Howev er , this decom- po sition is not appropriate for o ur purpos e b ecause the term 1234 stands for a non-simple co ne (pyramid), and hence this decomp osition does not simplify our pro blem. 5 Example 2.2. Consider the following system containing a redundant in- equality: + x 1 − x 2 + 0 x 3 ≤ 0 [1] − x 1 − x 2 + 0 x 3 ≤ 0 [2] 0 x 1 − x 2 + 0 x 3 ≤ 0 [3] The third inequalit y is redundan t. This example is ess en tially 2- dimensional, a nd 3 b oundaries (lines ) meet at the o rigin in R 2 . The simpli- cial co mplex is shown to b e F = { 1 , 2 , 3 , 12 , 13 , 23 , 1 23 } . This decomp osition is also unsuita ble beca use the 3-dimensional term 123 app ears althoug h this example is esse n tially 2-dimensional. T o av oid unfa vorable even ts such as have app eared in these e x amples we need the as sumption that { ∂ H i } is in the general p osition defined below. Definition 2.1. The set o f m hyper planes in R n , for example, { ∂ H i ⊂ R n | i = 1 , . . . , m } , is said to be in the genera l positio n when there do es not e x ist J ⊆ { 1 , . . . , m } s uc h that ∩ i ∈ J ∂ H i contains a max { n + 1 − | J | , 0 } dimensional affine subspace. Remark 2.1. This definition is equiv ale nt to that in Stanley (20 07). 7 The definition by Naiman a nd Wynn (1997) 6 is weak er than ours. In their defi- nition, when there do es not exist J ⊆ { 1 , . . . , m } , | J | = n + 1 , ∩ i ∈ J ∂ H i 6 = ∅ , the s et of hyperpla ne s is said to b e in the general p osition. Exa mple 2.2 is in the gener al p osition in their definition, wherea s it is not in the general po sition in our definition. When a given { ∂ H i } is not in the genera l p osition, we can apply an infinitesimal p erturbation to rearr ange the system into the gener al pos itio n. In this pap er we restrict our attention to the lexico graphic p erturbatio n in an outer direction prop osed b elow. F or b = ( b i ) ∈ R m and ε > 0 , define b ( ε ) = ( b i ( ε )) = ( b i + ε i ) = ( b 1 + ε 1 , . . . , b m + ε m ) ⊤ ( ε i is ε to the p ow er i ), and K ( ε ) = { x ∈ R n | A ⊤ x ≤ b ( ε ) } . (4) Define H i ( ε ), F i ( ε ) and F ( ε ) similarly . Note that K = K (0), H i = H i (0), F i = F i (0) and F = F (0 ). 6 The following tw o lemmas, Lemmas 2.1 a nd 2.2, are refinements o f Lemma 1 of Naiman and Wynn (1997) 6 under this lexicogra phic pe r turba- tion and with the stronger definition of the general p osition. Lemma 2.1. F or al l sufficiently smal l ε > 0 , (i) the family of hyp erplanes { ∂ H i ( ε ) | i = 1 , . . . , m } is in the gener al p osition, and (ii) F ( ε ) do es not dep end on ε . Write F (0+) = F ( ε ) . F (0+) is a simplicia l c omplex. Pro of. (i) Fix J ⊆ { 1 , . . . , m } . Consider first the case | J | ≤ n . Let A J = ( a i ) i ∈ J be n × | J | , and let b J ( ε ) = ( b i ( ε )) i ∈ J be | J | × 1 . Suppo se that ∩ i ∈ J ∂ H i ( ε ) contains an n +1 −| J | dimensiona l affine subspace { x = C d + x 0 | d ∈ R n +1 −| J | } , wher e C = ( c 1 , . . . , c n +1 −| J | ) is an n × ( n + 1 − | J | ) matrix with linearly indep enden t column v ectors. Then A ⊤ J ( C d + x 0 ) − b J ( ε ) = 0 for all d , and hence A ⊤ J C = 0 and A ⊤ J x 0 − b J ( ε ) = 0 hold. Since rank( C ) = n + 1 − | J | , we see tha t ra nk( A J ) ≤ | J | − 1. Hence A J can be written as A J = E F , where E is n × r and F is r × | J | with r = rank( A J ) ≤ | J | − 1. Substituting this, w e hav e 0 = A ⊤ J x 0 − b J ( ε ) = F ⊤ ( E ⊤ x 0 ) − b J ( ε ) = ( F ⊤ , − b J ( ε ))  E ⊤ x 0 1  . Therefore, f J ( ε ) := det  F − b J ( ε ) ⊤  ( F ⊤ , − b J ( ε ))  = det( F F ⊤ ) × b J ( ε ) ⊤  I − F ⊤ ( F F ⊤ ) − 1 F  b J ( ε ) = 0 . (5) Since rank( F ) = r , f J ( ε ) is a polynomia l in ε of at least degree 2 r ( ≥ 2). Then δ J = sup { δ > 0 | f J ( ε ) > 0 , ∀ ε ∈ (0 , δ ) } > 0 , since the n um be r of zeros of f J is finite. F or all ε ∈ (0 , δ J ), ∩ i ∈ J ∂ H i ( ε ) contains no ( n + 1 − | J | )-dimensional affine subspace. F or the case | J | ≥ n + 1 , as s ume that x 0 ∈ ∩ i ∈ J ∂ H i ( ε ) exists. Since r = rank( A J ) ≤ n ≤ | J | − 1, w e can follo w the pro of for | J | ≤ n ab ov e . W e define the po lynomial f J ( ε ) in (5), and conclude that for all ε ∈ (0 , δ J ), ∩ i ∈ J ∂ H i ( ε ) = ∅ . Let δ ∗ = min J δ J > 0 . F or all ε ∈ (0 , δ ∗ ), { ∂ H i ( ε ) | i = 1 , . . . , m } is in the g eneral p osition. 7 (ii) The pro o f for (i) implies that for a giv en J ⊆ { 1 , . . . , m } the fea si- bilit y or infeasibility is unchanged for a ll ε ∈ (0 , δ ∗ ). Corollary 2.1. F or e ach J ∈ F (0+ ) , (i) { a i | i ∈ J } is line arly indep endent, and (ii) the numb er of element s | J | is at most r ank( A ) . Pro of. (i) Note first that ∩ i ∈ J ∂ H i contains a p oint (i.e., 0-dimensio nal affine subspace), and becaus e { ∂ H i | i ∈ J } is in the general p osition (Lemma 2.1 , (i)), it should b e that ma x { n + 1 − | J | , 0 } ≥ 1 o r | J | ≤ n . Suppo se tha t { a i | i ∈ J } , J ∈ F (0+), is linearly dependent. Let L = span { a i | i ∈ J } and N = { x | a ⊤ i x = 0 , i ∈ J } . Then | J | ≥ dim( L ) + 1 = n − dim( N ) + 1 and hence dim( N ) ≥ n + 1 − | J | . On the other hand, ∩ i ∈ J ∂ H i ⊇ N + x 0 , w he r e x 0 ∈ ∩ i ∈ J ∂ H i . T his contradicts the fact that { ∂ H i | i ∈ J } is in the genera l p osition. (ii) The num b er of linearly indep enden t vectors { a i | i ∈ J } m ust b e at most dim span { a i } = rank( A ). Lemma 2. 2. L et F and F (0+) b e simplicial c omplexes define d in Pr op o- sition 2.1 and L emma 2. 1, r esp e ctively. Then F ⊇ F (0+) . Conse quent ly, |F | ≥ |F (0+) | . Pro of. Suppose that J ∈ F (0+). W e prov e that J ∈ F . F rom Lemma 2.1, for all s ufficien tly small ε > 0, it holds that E J ( ε ) ∩ K ( ε ) 6 = ∅ , where E J ( ε ) = ∩ i ∈ J ∂ H i ( ε ) = { x | A ⊤ J x = b ( ε ) } with A J = ( a j ) j ∈ J ( n × | J | ). Let E J = E J (0). F rom Corollar y 2.1, the column vectors of A J are linea rly independent, a nd hence E J 6 = ∅ . Assume tha t E J ∩ K = ∅ . Since b oth E J and K are clo sed p olyhedra in R n , there is a h ype r plane separa ting E J from K , a nd the distances of the hyperplane fr om E J and from K are gr eater than δ > 0 . By c ho osing ε > 0 such tha t max i =1 ,...,m ε k a i k < δ, this h yp erplane also separa tes E J ( ε ) from K ( ε ). There fo re, E J ( ε ) ∩ K ( ε ) = ∅ , which contradicts the a ssumption that J ∈ F (0+). Hence, E J ∩ K 6 = ∅ , from which J ∈ F follows. The following theorem gives another v ersion of the abstract tub e as- so ciated with a closed conv ex p olyhedron K . This is an improv emen t of 8 Corollar y 1 of Naiman and Wynn (1997), 6 who proved (6) for almos t a ll x (i.e., a weak abstra ct tub e). Theorem 2.1. The p air ( H , F (0+)) is an abstr act tub e in the sense that 1 ∪ m i =1 H c i ( x ) = X J ∈F (0+) ( − 1) | J |− 1 1 ∩ i ∈ J H c i ( x ) for al l x ∈ R n . (6) Pro of. As in Cor ollary 1 of Naiman a nd Wynn (1997), 6 it can b e ea sily prov ed that (6) holds for x / ∈ ∪ m i =1 ∂ H i . W e show that (6) holds for every x ∈ ∪ m i =1 ∂ H i . If x ∈ K ∩ ∪ m i =1 ∂ H i , i.e., x is on the boundar y of K , then x / ∈ H c i for all i , and (6) holds as 0 = 0. Suppo se that x ∈ K c ∩ ∪ m i =1 ∂ H i . Recall that for sufficiently small ε > 0, 1 ∪ m i =1 H i ( ε ) c ( x ) = X J ∈F (0+) ( − 1) | J |− 1 1 ∩ i ∈ J H i ( ε ) c ( x ) (7) holds. Since K is closed and ε is sufficient ly small, the left side of (7) bec omes 1 ∪ m i =1 H i ( ε ) c ( x ) = 1 ∪ m i =1 H c i ( x ) = 1. When x ∈ ∂ H i , a ⊤ i x = b i ≤ b i ( ε ), a nd hence 1 H c i ( x ) = 1 H i ( ε ) c ( x ) = 0. When x / ∈ ∂ H i , 1 H c i ( x ) = 1 H i ( ε ) c ( x ), becaus e the distanc e b et w een x and ∂ H i is pos itiv e. Therefore, in the right side of (7), 1 ∩ i ∈ J H i ( ε ) c ( x ) = Q i ∈ J 1 H i ( ε ) c ( x ) = Q i ∈ J 1 H c i ( x ) = 1 ∩ i ∈ J H c i ( x ). Summarizing the ab ov e, even when x ∈ K c ∩ ∪ m i =1 ∂ H i , the expr e ssion (7) with H i ( ε ) c replaced b y H c i holds, which is nothing but (6). Remark 2.2 . The a bstract tub e ( H , F (0 +)) substantially improv es the abstract tube ( H , F ) beca use (i) the num b er of elements of F (0+) is no t greater than that o f F (Lemma 2.2), and (ii) ∩ i ∈ J H c i , J ∈ F (0+ ), is a simple co ne in R n , or is the direct sum of a simple cone and a linear subspace (Co rollary 2.1). This is not alwa ys the c ase for F . W e present examples of abstr a ct tub es for p erturb ed polyhedr a b e lo w. Example 2.3. Applying a p erturbation to the ineq ua lit y sys tem in Ex am- 9 ple 2.1: − x 1 − x 2 + x 3 ≤ 1 + ε 1 [1] − x 1 + x 2 + x 3 ≤ 1 + ε 2 [2] + x 1 + x 2 + x 3 ≤ 1 + ε 3 [3] + x 1 − x 2 + x 3 ≤ 1 + ε 4 [4] where ε is an infinitesima l p ositive r eal num ber . The co r resp onding simpli- cial co mplex is shown to b e F (0+) = { 1 , 2 , 3 , 4 , 1 2 , 14 , 23 , 2 4 , 34 , 124 , 2 34 } , which is a prop er subset of F in (3). The num b er o f terms |F (0+) | = 11 is less than |F (0+) | = 15. The maximum num b er of elements of J ∈ F (0+) is 3 (= n ). The difference b etw een F and F (0 + ) is F \ F (0+) = { 13 , 123 , 134 , 123 4 } . The mem be r s of F \ F (0+) are c haracter ized as the sets containing 13 , that is, the non-existent edge in K ( ε ). Example 2.4. Applying a p erturbation to the inequality system in Exam- ple 2.2: + x 1 − x 2 + 0 x 3 ≤ 0 + ε 1 [1] − x 1 − x 2 + 0 x 3 ≤ 0 + ε 2 [2] 0 x 1 − x 2 + 0 x 3 ≤ 0 + ε 3 [3] The cor resp onding simplicial complex is shown to b e F (0+) = { 1 , 2 , 3 , 13 , 23 } . The lar gest | J | , J ∈ F (0 + ), is 2 , as pro v ed in (ii) of Corolla ry 2 .1 . Note that a different p erturbation obtained by alter ing the order o f inequalities + x 1 − x 2 + 0 x 3 ≤ 0 + ε 2 [1] − x 1 − x 2 + 0 x 3 ≤ 0 + ε 3 [2] 0 x 1 − x 2 + 0 x 3 ≤ 0 + ε 1 [3] yields a different simplicial complex F (0+) = { 1 , 2 , 12 } . This is an example wher e |F (0+) | dep ends on the p erturbation. 10 T aking the expe c ta tion of (6) with respect to the normally distr ibuted random v ector x , we obtain 1 − P ( K ) = X J ∈F (0+) ( − 1) | J |− 1 P ( ∩ i ∈ J H c i ) . As mentioned in Re ma rk 2.2, ∩ i ∈ J H c i app earing in the right side is of the class of cones whose probability can b e calculated by the metho d of Miw a et a l. (2003 ). 2 3. Construction of the abstract tub e In order to c onstruct the abstra ct tub e ( H , F (0+)), we ne e d to determine whether the s ystem a ⊤ i x = b i ( ε ) , i ∈ J, a ⊤ i x ≤ b i ( ε ) , i / ∈ J, has a so lution for a given subset J ⊆ { 1 , . . . , m } , whe r e ε is an infinitesima l po sitiv e real n um ber . F o r that purp ose, linear progra mming (LP) techniques for deg enerate sy stems are useful. First r ewrite the problem in a linear progra mming forma t as follows: Problem 3.1. F or a given subset J ⊆ { 1 , . . . , m } determine whether ther e exists a fe asible s olut ion x, y ≥ 0 ( x, y ∈ R n ) and u i ≥ 0 ( i = 1 , . . . , m ) , v j ≥ 0 ( j ∈ J ) su ch that u i + a ⊤ i ( x − y ) − b i ( ε ) = 0 , i = 1 , . . . , m, v j − a ⊤ j ( x − y ) + b j ( ε ) = 0 , j ∈ J. T o solve this problem w e first define a tableau (matr ix ) corr espo nding to the a b ov e linear system. Let A J = ( a j ) j ∈ J ( n × | J | ), b J ( ε ) = ( b j ( ε )) j ∈ J ( | J | × 1), and define a n M × N ma trix [ t ij ] M × N =   A ⊤ − A ⊤ − b ( ε ) − A ⊤ J A ⊤ J b J ( ε ) 1 ⊤ n 1 ⊤ n 0   M × N where M = m + | J | + 1, N = 2 n + 1 a nd 1 ⊤ n = (1 , . . . , 1) 1 × n . An algorithm for checking the feasibility of J is giv en as fo llo ws (Theorem 4-4-3 of Ir i, 1973 8 ). Algorithm 3.1. 11 1. Cho ose an index i (=: i 0 ) such that i ≤ M − 1 and t iN > 0 . If no such i exists, then the system is fe asible. 2. Among j ’s su ch that j ≤ N − 1 and t i 0 ,j < 0 (if no such j ex- ists, t hen the system is not fe asible), cho ose an index j (=: j 0 ) t hat maximizes t M j / | t i 0 ,j | . 3. F or i ∈ { 1 , . . . , M } \ { i 0 } and j ∈ { 1 , . . . , N } \ { j 0 } , let t ij := t ij − t ij 0 t i 0 j /t i 0 j 0 , t ij 0 := t ij 0 /t i 0 j 0 , t i 0 j := − t i 0 j /t i 0 j 0 , t i 0 j 0 := 1 /t i 0 j 0 . 4. Go t o Step 1. In Step 1, t iN is a p olynomial in ε . t iN > 0 me ans that k (0 ≤ k ≤ n ) exists such that t iN = c k ε k + c k +1 ε k +1 + · · · + c n ε n and c k > 0 . Using Algorithm 3.1, for every nonempty subset J ⊆ { 1 , . . . , m } we can confirm whether J is feasible (i.e., J ∈ F (0 +)). The following p oints should be incor por ated int o the implementation of the a lgorithm. Remark 3.1. Thanks to Corollar y 2 .1, the ra ng e of J for chec king the feasibility can b e restricted to J r,m = { J | J ⊆ { 1 , . . . , m } , 0 < | J | ≤ r } , (8) where r = rank( A ). Mo reov e r , if J 1 ⊂ J 2 and J 1 / ∈ F (0 +), then J 2 / ∈ F (0 +). This re duce s the range for checking the feasibility . Remark 3. 2 . In Algo r ithm 3 .1 the p o lynomial t iN = c 0 + c 1 ε 1 + · · · + c n ε n can b e repres e n ted a s an n + 1 vector ( c 0 , c 1 , . . . , c n ). In this repre s en tation t iN > 0 means ( c 0 , c 1 , . . . , c n ) > (0 , . . . , 0) in the lexicog r aphic or der. T his techn ique is known as the lexicogra phic metho d in linea r prog ramming. 8,9 Remark 3.3. In Algorithm 3.1 the judgment that a v alue z is p ositive ( z > 0) should be implemen ted as z > eps , where eps is a small po sitive n um be r comparable to the machine epsilon. F or example, the ma c hine epsilon in R (32-bit) is 2 − 52 . = 2 × 10 − 16 , and in our proto t ype R progr am eps is set to 10 − 14 . 12 Remark 3.4 . In c o ding a computer pro gram it is co n venien t to treat the candidate s et J r,m in (8) as a totally or dered set. Note that |J r,m | = r X j =1  m j  , |J 0 ,m | = 0 . W e intro duce a total order as follows: F or J 1 = { d 1 , . . . , d l 1 } ( d 1 < · · · < d l 1 ) a nd J 2 = { e 1 , . . . , e l 2 } ( e 1 < · · · < e l 2 ) in J r,m , define J 1  J 2 if l 1 < l 2 , or l 1 = l 2 and ( d 1 , . . . , d l 1 ) ≤ ( e 1 , . . . , e l 1 ) (lexicogra phically). Then J = { d 1 , . . . , d l } ( d 1 < · · · < d l ) is the L -th smallest element in J r,m , where L = 1 + |J l − 1 ,m | + l X j =1 d j − 1 X k = d j − 1 +1  m − k l − j  . Here we le t d 0 = 0, a nd P d j − 1 k = d j − 1 +1 = 0 if d j − 1 + 1 > d j − 1. Conv er sely , for a given L suc h that 1 ≤ L ≤ |J r,m | , the corre s ponding L -th smallest element J = { d 1 , . . . , d l } ( d 1 < · · · < d l ) in J r,m is obtained as fo llows: l := max n l ≥ 1   |J l − 1 ,m | < L o , and for j = 1 , . . . , l , itera tiv ely , d j := max  d j ≥ d j − 1 + 1    j X i =1 d i − 1 X k = d i − 1 +1  m − k l − i  < L − |J l − 1 ,m |  . 4. T uk ey’s studen tized range The multiple c o mparisons pr oc e dure is a typical a pplication where the mul- tidimensional normal probability in (1) is req uir ed. This is a sta tistical pro - cedure for com bining several statistica l tests. As an exa mple, consider the compariso n of k normal means. Suppose that, for i = 1 , . . . , k , X i is inde- pendently distributed according to a normal distribution with mea n µ i and v ariance σ 2 i , a nd consider testing the equality of the means µ i , i = 1 , . . . , k . F or X = ( X 1 , . . . , X k ), T ukey’s studentized r ange sta tistic is defined as T ( X ) = max 1 ≤ i 0 . O nc e an abstra ct tub e is obtained for some c > 0, this abstra ct tub e is av a ilable for computing the distribution function P r( T ( X ) ≤ c ) for any c in the n ull a nd non-null cases . −4 −2 0 2 4 0.950 0.960 0.970 s prob Fig. 2. T ukey-Kramer conject ure ( k = 5) References 1. A. Barvinok, A Course in Convexity (AMS , 2002). 2. T. Miwa, A. J. H a yter and S. Ku riki, J. R. Statist. So c., B 65 , 223 (2003). 3. I. G. Ab rahamson, Ann. Math. Statist. 35 , 1685 (1964). 4. A. Genz, F. Bretz, T. Miw a, X. Mi, F. Leisch, F. Scheipl and T. Hothorn, mvtnorm: Mul tivariate Normal and t D istributions , (2011). R pack age version 0.9-9991. 5. A. Genz and F. Bretz, Computation of Multivariate Normal and t Pr ob abili- ties, L e ctur e Notes in Statistics (Springer, Heidelberg, 2009). 6. D. Q . Naiman and H. P . Wynn, A nn. Statist. 25 , 1954 (1997). 7. R . P . Stanley , An intr o duction to hyp erplane arr angements , in Ge ometric 15 c ombinatorics , eds. V. R. Ezra Miller and B. S turmfels, IAS/Park Cit y Math- ematics, V ol. 13 (American Mathematical So ciety , 2007), pp. 389–496. 8. M . Iri, Line ar Pr o gr amming (in Jap anese) (Hakujitsusha, 1973). 9. R. J. V anderb ei, Line ar Pr o gr amming: F oundations and Extensions , second edn. (Kluw er, 2001). 10. A. J. Ha yter, Ann. Statist. 12 , 61 (1984).