Holonomic gradient method for the probability content of a simplex region with a multivariate normal distribution
aa r X i v : . [ m a t h . S T ] D ec Holonomic gradient method for the probabilitycontent of a simplex region with a multivariatenormal distribution
Tamio Koyama
Abstract
We use the holonomic gradient method to evaluate the probabilitycontent of a simplex region with a multivariate normal distribution. Thisprobability equals to the integral of the probability density function ofthe multivariate Gaussian distribution on the simplex region. For thispurpose, we generalize the inclusion–exclusion identity which was givenfor polyhedra, to the faces of a polyhedron. This extended inclusion–exclusion identity enables us to calculate the derivatives of the functionassociated with the probability content of a polyhedron in general position.we show that these derivatives can be written as integrals of the faces ofthe polyhedron.
The holonomic gradient method (HGM) is an algorithm for the numerical cal-culation of holonomic functions. It is a variation on the holonomic gradientdescent (HGD) proposed in [8]. A holonomic function is an analytic function ofseveral variables which satisfies a holonomic system. Here, a holonomic systemrefers to a system of linear differential equations with polynomial coefficientswhich induces a holonomic module in terms of D -module theory [9]. The HGMevaluates a holonomic function by numerically solving an initial value problemfor an ordinary differential equation. This ordinary differential equation is de-rived from the Pfaffian equation (an integrable connection) associated with thefunction. For details, see [3] and its references. Several normalizing constantsand the probability content of a region can be regarded as a holonomic functionwith respect to their parameters, and we can use the HGM to estimate thesolution to this function. For example, the HGM was used to evaluate the cu-mulative distribution function for the largest root of the Wishart matrix in [2],and we utilized the HGM for the orthant probability and the Fisher–Binghamintegral in [4] and [5], respectively.Our motivation is to apply the HGM to the numerical calculation of theprobability content of a polyhedron with a multivariate normal distribution. Apolyhedron is a subset of a d -dimensional Euclidean space R d which is defined by1 finite number of linear inequalities. Intervals of real numbers, orthants, cubes,and simplexes are examples of polyhedra. In [7], Naiman and Wynn describedexamples of how the evaluation of the probability content of a polyhedron canbe used to find critical probabilities for multiple comparisons.The probability content of a polyhedron is a generalization of orthant prob-abilities discussed in [4] since the orthant probabilities can be expressed as theprobability content of a simplicial cone. We derived the HGM for orthant prob-abilities in [4], and its implementation in [6]. A study of phylogenomics utilizedour implementation in [11].In order to utilize the HGM for the numerical calculation of the probabilitycontent of a polyhedron, we need to provide the Pfaffian equation explicitly andto evaluate the initial value for the ordinary differential equation. Our previouspaper [3] showed that the probability content of a polyhedron in general positioncan be expressed as an analytic function and we explicitly provided a holonomicsystem and Pfaffian equations for this function. In this paper, we calculatethe derivatives of the function, and show that these derivatives can be writtenas integrals of the faces of the polyhedron. This result provides formulae tocompute the initial value exactly for the cases where the polyhedron is in generalposition and bounded, or the polyhedron is a simplicial cone.In order to calculate the derivatives, we generalize the inclusion–exclusionidentity that was given for polyhedra in [1], to the faces of the polyhedron.Since a face of a polyhedron is also a polyhedron, the inclusion–exclusion iden-tity for the face holds, i.e., the indicator function of the face can be writtenas a linear combination of indicator functions of simplicial cones. Our general-ized inclusion–exclusion identity gives the expression for this linear combinationexplicitly.In the numerical experiments, we consider simplexes and simplicial cones, asthese are fundamental examples of polyhedra. Utilizing the theoretical resultsconcerned with the Pffafian equation and the initial value, we implement theHGM to evaluate the probability contents of a simplex and a simplicial cone.We show that our implementation works well for a 10-dimensional simplex.This paper is organized as follows. In section 2 we review results fromour previous paper [3]. In section 3 we extend the inclusion–exclusion identitywhich was given for polyhedra in [1], and provide an analogous formula forthe indicator function of a face of a polyhedron. In section 4 we calculate thederivatives of the function defined by the probability content of a polyhedronfor the multivariate normal distribution, and show that these derivatives can bewritten by integrals on corresponding faces. In section 5 attention is directedtowards the case where the polyhedron in general position is bounded, and thecase where the polyhedron is a simplicial cone. We discuss the evaluation of themultivariate normal probabilities of polyhedra by the HGM in these two cases.In section 6 we present numerical examples. Acknowledgements.
This work was supported by JSPS KAKENHI GrantNumber 263125. 2
Summary of Previous Work
In this section we review the results of our previous paper [3]. Let us considera polyhedron P := ( x ∈ R d : d X i =1 ˜ a ij x i + ˜ b j ≥ , ≤ j ≤ n ) (1)where ˜ a ij , ˜ b j (1 ≤ i ≤ d, ≤ j ≤ n ) are real numbers. We denote by ˜ a and ˜ b the d × n matrix (˜ a ij ) and the vector (˜ b , . . . , ˜ b n ) ⊤ respectively. We suppose thatthe polyhedron P is in general position and its bounding half-spaces are H j := ( x ∈ R d : d X i =1 ˜ a ij x i + ˜ b j ≥ ) (1 ≤ j ≤ n ) . (2)For the definitions of general position and the set of the bounding half-spacesfor a polyhedron, see [3].We denote by F the abstract simplicial complex associated with the poly-hedron P , i.e., F := J ⊂ { , , . . . , n } | \ j ∈ J H j = ∅ . Let ϕ ( a, b ) = Z R d π ) d/ exp − d X i =1 x i ! X F ∈F Y j ∈ F ( H ( f j ( a, b, x )) − dx (3)be a function with variables a ij , b j (1 ≤ i ≤ d, ≤ j ≤ n ). Here, H ( x ) is theHeaviside function and we set f j ( a, b, x ) = P di =1 a ij x i + b j . We denote by a and a j ( j = 1 , dotsn ) the d × n matrix with elements a ij , ( i = 1 , . . . d ) and the columnvector ( a j , . . . , a dj ) ⊤ respectively. And b is a column vector ( b , . . . , b n ) ⊤ withlength n . For J ∈ F , we put g J ( a, b ) = Y j ∈ J ∂ b j • ϕ ( a, b ) . (4)And let g ( a, b ) = ( g J ( a, b )) J ∈F be a vector-valued function, then g ( a, b ) satisfiesthe following Pfaffian equations [3, Theorem 22] : ∂ a ij g J = n X k =1 a ik ∂ b k ∂ b j g J (1 ≤ i ≤ d, ≤ j ≤ n, J ∈ F ) , (5) ∂ b j g J = g J ∪{ j } ( j ∈ J c , J ∈ F ) , (6) ∂ b j g J = − X k ∈ J α jkJ ( a ) b k g J + X ℓ ∈ J c α kℓ ( a ) g J ∪ ℓ ! ( j ∈ J, J ∈ F ) . (7)3ere, ( α ijF ( a )) i,j ∈ F is the inverse matrix of α F ( a ) = (cid:16)P dk =1 a ki a kj (cid:17) i,j ∈ F , whichis a submatrix of the gram matrix of a . Note that the right hand side of (5) canbe rewritten, with recourse to (6) and (7), as a linear combination of g J withrational functions as coefficients. Let P be the polyhedron defined by (1), and suppose the family of the boundinghalf-spaces for P is given by (2). Then, we have the following inclusion–exclusionidentity [1]: n Y j =1 H d X i =1 ˜ a ij x i + ˜ b j ! = X J ∈F Y j ∈ J H d X i =1 ˜ a ij x i + ˜ b j ! − ! ( x ∈ R d ) . (8)In [3], we showed that if the polyhedron P is in general position, there exists aneighborhood U of the parameter (˜ a, ˜ b ) ∈ R d × n × R n such that n Y j =1 H d X i =1 a ij x i + b j ! = X J ∈F Y j ∈ J H d X i =1 a ij x i + b j ! − ! (9)holds for any ( a, b ) ∈ U and x ∈ R d . The left hand sides of (8) and (9) arethe indicator functions for the corresponding polyhedra. In this section we giveanalogous identities for the indicator functions of a face of the polyhedra.Let F be the abstract simplicial complex for the polyhedron P . For J ∈ F , F J := { F ∈ F | J ⊂ F } . For parameter ( a, b ) ∈ R d × n × R n and J ∈ F , we define a hyperplane V ( J, a, b )by V ( J, a, b ) = ( x ∈ R d | d X i =1 a ij x i + b j = 0 ( j ∈ J ) ) . (10) Proposition 1.
Suppose the polyhedron P is in general position. For each J ∈ F , we have the equation Y j ∈ [ n ] \ J H d X i =1 ˜ a ij x i + ˜ b j ! = X F ∈F J Y j ∈ F \ J H d X i =1 ˜ a ij x i + ˜ b j ! − ! . for any x ∈ V ( J, ˜ a, ˜ b ) .Proof. Let s be the number of elements in the set J . Since the polyhedron P is in general position, we have s ≤ d . By replacing the indices, we can assume J = { n − s + 1 , . . . , n } without loss of generality. Applying the Euclideantransformation for P , we can assume ˜ a ij = 0 , ˜ b j = 0 (1 ≤ i ≤ d − s, n − s + 1 ≤ ≤ n ). Then, by the assumption of general position, the vectors ˜ a n − s +1 , . . . , ˜ a n are linearly independent (see, [3, Corollary 20]). Hence we have V ( J, ˜ a, ˜ b ) = { x ∈ R d | x d − s +1 = · · · = x d = 0 } , and the problem is reduced to the proof of n − s Y j =1 H d − s X i =1 ˜ a ij y i + ˜ b j ! = X F ∈F J Y j ∈ F \ J H d − s X i =1 ˜ a ij y i + ˜ b j ! − ! (11)for arbitrary y = ( y , . . . , y d − s ) ⊤ ∈ R d − s . Let us consider a polyhedron P ′ := ( y ∈ R d − s | d − s X i =1 ˜ a ij y i + ˜ b j ≥ , ≤ j ≤ n − s ) . Suppose that there are t redundant inequalities in the definition of P ′ . By replac-ing indices, we can assume that the redundant inequalities are P d − si =1 ˜ a ij y i +˜ b j ≥ , ( n − s − t + 1 ≤ j ≤ n − s ) . Then, all facets of P ′ are given by F ′ j := P ′ ∩ ( y ∈ R d − s | d − s X i =1 ˜ a ij y i + ˜ b j = 0 ) (1 ≤ j ≤ n − s − t ) , and the abstract simplicial complex for P ′ is F ′ = J ′ ⊂ { , , . . . , n − s − t } | \ j ∈ J ′ F ′ j = ∅ . Applying the inclusion-exclusion identity for P ′ , we have n − s − t Y j =1 H d − s X i =1 ˜ a ij y i + ˜ b j ! = X F ∈F ′ Y j ∈ F H d − s X i =1 ˜ a ij y i + ˜ b j ! − ! . (12)Since the left hand sides of (11) and (12) are both equal to the indicatorfunction of P ′ , they are equal to each other. Consequently, we need to showthat the right hand sides of (11) and (12) are equal. It is easy to show that themapping ψ : F ′ → F J ( J ′ J ∪ J ′ )is a bijection. Rewriting the right hand side of (12) in terms of F J , we have thesame expression for the right hand side of (11).For use in the next section, we extend Proposition 1. We introduce the5ollowing notation. For parameters ( a, b ) ∈ R d × n × R n , let H j ( a, b ) = ( x ∈ R d | d X i =1 a ij x i + b j ≥ ) (1 ≤ j ≤ n ) , H ( a, b ) = {H ( a, b ) , . . . , H n ( a, b ) } ,P ( a, b ) = n \ j =1 H j ( a, b ) ,F j ( a, b ) = ( x ∈ R d | d X i =1 a ij x i + b j = 0 ) ∩ P ( a, b ) , F ( a, b ) = J ⊂ [ n ] | \ j ∈ J F j ( a, b ) = ∅ . The following lemma holds.
Lemma 1.
Suppose the polyhedron P is in general position. Then, there existsa neighborhood U of the parameter (˜ a, ˜ b ) ∈ R d × n × R n which satisfies the follow-ing: for any parameter ( a, b ) in U , P ( a, b ) is in general position and F ( a, b ) = F holds.Proof. We put a i = 0 (1 ≤ i ≤ d ), b = 0, andˆ F j ( a, b ) = ( x , x ) ∈ R × R d | P di =1 a ij x i + b j = 0 , P di =1 a ik x i + b k ≥ ≤ k ≤ n ) (0 ≤ j ≤ n )ˆ F ( a, b ) = J ⊂ { , , . . . , n } | \ j ∈ J ˆ F j ( a, b ) = ∅ . By Theorem 23 in [3], the set U := (cid:26) ( a, b ) ∈ R d × n × R n | P ( a, b ) is in general position , ˆ F ( a, b ) = ˆ F (˜ a, ˜ b ) (cid:27) is a neighborhood of the point (˜ a, ˜ b ) ∈ R d × n × R n . Consider arbitrary ( a, b ) ∈ U ,from Corollary 19 in [3], we have F ( a, b ) = F (˜ a, ˜ b ) = F from ˆ F ( a, b ) = ˆ F (˜ a, ˜ b ).By Lemma 22 in [3], all facets of P ( a, b ) are given by F j ( a, b ) ( { j } ∈ F ( a, b )).The equation F ( a, b ) = F implies { j } ∈ F ( a, b ) for all j ∈ [ n ]. Consequently, H ( a, b ) is the bounding half-spaces for P ( a, b ) and P ( a, b ) is in general position.Finally, we have the following. 6 heorem 1. Suppose the polyhedron P is in general position and J ∈ F . Thereexists a neighborhood U of (˜ a, ˜ b ) ∈ R d × n × R n such that for any ( a, b ) ∈ U and x ∈ V ( J, a, b ) , we have Y j ∈ [ n ] \ J H d X i =1 a ij x i + b j ! = X F ∈F J Y j ∈ F \ J H d X i =1 a ij x i + b j ! − ! . (13) Proof.
Let U be a neighborhood of (˜ a, ˜ b ) in Lemma 1. Then the polyhedron P ( a, b ) is in general position. By Lemma 22 in [3], the abstract simplicialcomplex associated with P ( a, b ) is equivalent to F ( a, b ). The equation F ( a, b ) = F implies J ∈ F ( a, b ). Hence, we can apply Proposition 1 which gives Y j ∈ [ n ] \ J H d X i =1 a ij x i + b j ! = X F ∈F J ( a,b ) Y j ∈ F \ J H d X i =1 a ij x i + b j ! − ! for any x ∈ V ( J, a, b ). Here, we put F J ( a, b ) = { F ⊂ F ( a, b ) | J ⊂ F } . Since F ( a, b ) = F implies F J ( a, b ) = F J , we thus have equation (13). In this section we derive an expression for the function g J ( a, b ) ( J ∈ F ), whichis a derivative of the function ϕ ( a, b ) defined by the probability content of thepolyhedron P for the multivariate normal distribution. We then show that thefunction g J ( a, b ) can be expressed as an integral on the hyperplane (10). Forsimplicity, we put ϕ F ( a, b ) = Z R d exp − d X i =1 x i ! Y j ∈ F H ( − f j ( a, b, x )) dx. (14)Then, the function ϕ ( a, b ) in (3) can be written as ϕ ( a, b ) = X F ∈F ( − | F | (2 π ) d/ ϕ F ( a, b ) . In order to obtain expressions for the function g J ( a, b ), we first consider ∂ Jb • ϕ F ( a, b ) for F ∈ F . Consider the following case: Lemma 2.
Let ≤ p ≤ q ≤ d , J = { , . . . , p } , and F = { , . . . , p, . . . , q } .Suppose that parameter a satisfies a ij = 0 ( p < i ≤ d, ≤ j ≤ p ) and α F ( a ) =7 P dk =1 a ki a kj (cid:17) i,j ∈ F is a non-singular matrix, i.e., a = ( a ij ) = a · · · a p a p +1) · · · a q ∗ · · · ∗ ... ... ... ... ... ... a p · · · a pp a p ( p +1) · · · a pq ∗ · · · ∗ · · · a ( p +1)( p +1) · · · a ( p +1) q ∗ · · · ∗ ... ... ... ... · · · a d ( p +1) · · · a dq ∗ · · · ∗ and the vectors a , . . . , a p , . . . , a q are linearly independent. Then, the function ∂ Jb • ϕ F ( a, b ) is equal to the integral ( − | J | p | α J ( a ) | Z V ( J,a,b ) exp − d X i =1 x i ! Y j ∈ F \ J H ( − f j ( a, b, x )) µ ( dx ) . (15) Here, µ is the volume element of the hyperplane V ( J, a, b ) .Proof. Let a d × d matrix U = ( u ij ) be U = a · · · a p · · · a p · · · a pp · · · · · · · · · . The elements of U can be written as u ij = ( a ij (1 ≤ j ≤ p ) δ ij ( p < j ≤ d ) . Here, δ ij is Kronecker’s delta. As the vectors a , . . . , a p are linearly independent,the matrix U is regular, and we have | U | = | α J ( a ) | . We denote the inversematrix of U by U − = ( u ij ). Consider a transformation of variables y j = P di =1 u ij x i (1 ≤ j ≤ d ) . By the relationships x i = (P pk =1 u ki y k (1 ≤ i ≤ p ) y i ( p + 1 ≤ i ≤ d ) ,y j = d X i =1 a ij x i + b j (1 ≤ j ≤ p ) , ϕ F ( a, b ) can be written as1 p | α J ( a ) | Z R d e − P pi =1 ( P pk =1 u ki y k ) − P di = p +1 y i p Y j =1 H ( − y j − b j ) × q Y j = p +1 H − p X i =1 p X k =1 a ij u ki y k − d X i = p +1 a ij y i − b j dy . . . dy d = Z − b −∞ · · · Z − b p −∞ G ( y , . . . , y p ; a, b, U ) dy p . . . dy . Here, we put G ( y , . . . , y p ; a, b, U )= e − P pi =1 ( P pk =1 u ki y k ) p | α J ( a ) | Z ∞−∞ · · · Z ∞−∞ e − P di = p +1 y i × q Y j = p +1 H − p X i =1 p X k =1 a ij u ki y k − d X i = p +1 a ij y i − b j dy p +1 . . . dy d . Since G ( y , . . . , y p ; a, b, U ) is a continuous function with respect to y , . . . , y p ,we have ∂ Jb • ϕ F ( a, b ) = ( − p G ( − b , . . . , − b p ; a, b, U ).When p = d , we have ∂ Jb • ϕ F ( a, b ) = ( − d p | α J ( a ) | e − P di =1 ( P dk =1 u ki b k ) Since − ( U − ) ⊤ b is the unique point in V ( J, a, b ), ∂ Jb • ϕ F ( a, b ) equals to (15).Suppose p = d , and define a mapping ψ ( x ) = ( y p +1 , . . . , y d ) for the hyper-plane V ( J, a, b ) to R d − p by y j = x j ( p + 1 ≤ j ≤ d ), then this mapping is a localcoordinate system on V ( J, a, b ). At this coordinate, the functions x , . . . , x d on V ( J, a, b ) can be written as x i = ( − P pk =1 u ki b k (1 ≤ i ≤ p ) ,y i ( p + 1 ≤ i ≤ d ) . The Riemannian metric induced on V ( J, a, b ) is P di =1 dx i ⊗ dx i = P dj = p +1 dy j ⊗ dy j . Calculating (15) with this coordinate, we have g J ( a, b ) = ( − p G ( − b , . . . , − b p ; a, b, U ).We now extend this lemma. Lemma 3.
Let ≤ p ≤ q ≤ d , J = { , . . . , p } , and F = { , . . . , p, . . . , q } .Suppose α F ( a ) is a regular matrix. Then, the function ∂ Jb • ϕ F ( a, b ) is equal to (15) . roof. It is sufficient to show that this reduces to Lemma 2.For a suitable special orthogonal matrix R , the d × n matrix a ′ := Ra satisfiesthe condition a ′ ij = 0 ( p < i ≤ d, ≤ j ≤ p ). Since α F ( a ) = α F ( a ′ ), α F ( a ′ ) isalso a regular matrix by this assumption. Hence, the parameter ( a ′ , b ) satisfiesthe condition of Lemma 2.Since the Lebesgue measure is invariant under the action of the special or-thogonal group, we have ϕ F ( a, b ) = ϕ F ( a ′ , b ) for any b ∈ R n . Consequently, wehave ∂ Jb • ϕ F ( a, b ) = ∂ Jb • ϕ F ( a ′ , b ), and 1 / p | α J ( a ) | = 1 / p | α J ( a ′ ) | .Considering (15), we put˜ ϕ F ( a, b ) = Z V ( J,a,b ) exp − d X i =1 x i ! Y j ∈ F \ J H ( − f j ( a, b, x )) µ ( dx ) . We need to show ˜ ϕ F ( a, b ) = ˜ ϕ F ( a ′ , b ). When p = d , this relation is trivial.Suppose that p < d . Take vectors v j = ( v j , . . . , v dj ) ⊤ ( q + 1 ≤ j ≤ d ) suchthat a , . . . , a q , v q +1 , . . . , v d are linearly independent. Let U = ( u ij ) be a matrixobtained by arranging these vectors, u ij = ( a ij (1 ≤ j ≤ q ) ,v ij ( q + 1 ≤ j ≤ d ) . We denote the inverse matrix of U , by U − = ( u ij ). And define a matrix U ′ = ( u ′ ij ) as U ′ = RU , and denote its inverse by U ′− = ( u ′ ij ).First, we calculate ˜ ϕ F ( a, b ). If we define a map ψ ( x ) = ( y p +1 , . . . , y d ) fromthe hyperplane V ( J, a, b ) to R d − p by y j = P di =1 u ij x i ( p + 1 ≤ j ≤ d ), then itis a local coordinate system on V ( J, a, b ). With this coordinate, the function x i on V ( J, a, b ) can be written as x i = − p X k =1 u ki b k + d X k = p +1 u ki y k (1 ≤ i ≤ d ) . Hence, the Riemannian metric on the hyperplane V ( J, a, b ) is d X i =1 dx i ⊗ dx i = p X k =1 p X ℓ =1 d X i =1 u ki u ℓi ! dy k ⊗ dy ℓ . Let D be the determinant of the matrix (cid:16)P di =1 u ki u ℓi (cid:17) ≤ k,ℓ ≤ p . The integral˜ ϕ F ( a, b ) can be written as1 p | D | Z R d − p e − P di =1 ( − P pk =1 u ki b k + P dk = p +1 u ki y k ) q Y j = p +1 H ( − y j − b j ) d Y j = p +1 dy j . Next, we calculate ˜ ϕ F ( a ′ , b ). If we define a map ψ ′ ( x ) = ( y p +1 , . . . , y d ) fromthe hyperplane V ( J, a ′ , b ) to R d − p by y j = P di =1 u ′ ij x i ( p + 1 ≤ j ≤ d ), then it10s a local coordinate system on V ( J, a ′ , b ). With this coordinate, the function x i on V ( J, a ′ , b ) can be written as x i = − p X k =1 u ′ ki b k + d X k = p +1 u ′ ki y k (1 ≤ i ≤ d ) . By U ′− = U − R ⊤ , the Riemannian metric on V ( J, a ′ , b ) is d X i =1 dx i ⊗ dx i = p X k =1 p X ℓ =1 d X i =1 u ′ ki u ′ ℓi ! dy k ⊗ dy ℓ = p X k =1 p X ℓ =1 d X i =1 u ki u ℓi ! dy k ⊗ dy ℓ . And we have d X i =1 x i = d X i =1 − p X k =1 u ′ ki b k + d X k = p +1 u ′ ki y k = d X i =1 − p X k =1 u ki b k + d X k = p +1 u ki y k . Hence we have ˜ ϕ F ( a ′ , b ) = ˜ ϕ F ( a, b ).From Lemma 3, we have the following. Lemma 4.
Let F ∈ F and suppose α F ( a ) is a regular matrix. Then, we have ∂ Jb • ϕ F ( a, b )= ( − | J | √ | α J ( a ) | R V ( J,a,b ) e − P di =1 x i Q j ∈ F \ J H ( − f j ( a, b, x )) µ ( dx ) ( J ⊂ F ) , J F ) . Proof.
When J ⊂ F , it reduces to Lemma 3 since we can assume J = { , . . . , p } , F = { , . . . , p, . . . , q } , and 1 ≤ p ≤ q ≤ d without loss of generality. When J F , the integral in (14) does not depend on the variables b j ( j ∈ J \ F ).Hence, the derivative with respect to b j is 0. Theorem 2.
Suppose that the polyhedron P is in general position and J ∈ F ,then there exists a neighborhood U of the parameter (˜ a, ˜ b ) ∈ R d × n × R n suchthat the equation g J ( a, b ) = 1(2 π ) d/ p | α J ( a ) | Z V ( J,a,b ) e − P di =1 x i Y j ∈ [ n ] \ J H ( f j ( a, b, x )) µ ( dx )(16) holds for any ( a, b ) ∈ U . roof. By (3) and (4), we have(2 π ) d/ g J ( a, b ) = X F ∈F ∂ Jb • ( − | F | Z R d e − P di =1 x i Y j ∈ F H ( − f j ( a, b, x )) dx. Applying Lemma 4 to each term on the right hand side of the above equation,we can show that (2 π ) d/ g J ( a, b ) is equal to1 p | α J ( a ) | Z V ( J,a,b ) e − P di =1 x i X F ∈F J ( − | F \ J | Y j ∈ F \ J H ( − f j ( a, b, x )) µ ( dx )= 1 p | α J ( a ) | Z V ( J,a,b ) e − P di =1 x i X F ∈F J Y j ∈ F \ J H ( f j ( a, b, x ) − µ ( dx ) . From Theorem 1, we have the equation (16).
In this section, we discuss the computation of the probability content of a poly-hedron with a multivariate normal distribution for the case where the polyhe-dron is in general position and bounded, and the case where the polyhedron isa simplicial cone.
Let us consider the case where the polyhedron P in general position is bounded. Lemma 5.
Suppose the polyhedron P is bounded. Then, the set ( x ∈ R d | d X i =1 ˜ a ij x i ≥ ≤ j ≤ n ) ) (17) contain only the origin.Proof. By Proposition 1.12 in [12], the set (17) is equal to (cid:8) y ∈ R d | x + ty ∈ P ( x ∈ P, t ≥ (cid:9) . Since P is bounded, this set does not contain any element except the origin. Proposition 2.
Suppose the polyhedron P in general position is bounded. Then,for J ∈ F , we have g J (˜ a,
0) = ( √ | α J (˜ a ) | ( | J | = d )0 ( else ) . (18)12 roof. Calculating the left hand side, we have(2 π ) d/ g J (˜ a, t → +0 (2 π ) d/ g J (˜ a, tb )= lim t → +0 p | α J (˜ a ) | Z V ( J, ˜ a,tb ) e − P di =1 x i Y j ∈ [ n ] \ J H ( f j (˜ a, tb, x )) µ ( dx )= 1 p | α J (˜ a ) | Z V ( J, ˜ a, e − P di =1 x i Y j ∈ [ n ] \ J H ( f j (˜ a, , x )) µ ( dx ) . By Lemma 5, the integral domain is { } . Hence we have (18).Consequently, in order to compute the probability content of P for a multi-variate normal distribution, we can take the path of the HGM as a ( t ) = ˜ a, b ( t ) = ˜ tb (0 ≤ t ≤ . This path does not pass through the singular locus of the Pfaffian equations (5),(6) and (7). The initial value g J ( a (0) ,
0) at t = 0 is given explicitly by (18). Consider the case where the polyhedron P is a simplicial cone, i.e., n = d andthe vectors ˜ a , . . . , ˜ a d are linearly independent. We can assume without loss ofgenerality that ˜ a is an upper triangular matrix. Then define γ ( t ) = ( a ( t ) , b ( t ))by a ( t ) = (1 − t )diag(˜ a , . . . , ˜ a dd ) + t ˜ a, b ( t ) = t ˜ b (0 ≤ t ≤ . This does not pass through the singular locus of the Pfaffian equation. Theinitial value is g J ( a (0) , b (0)) = 1 (cid:12)(cid:12)(cid:12)Q j ∈ J ˜ a jj (cid:12)(cid:12)(cid:12) r π d −| J | . In this section we compare the performance of our HGM method with a MonteCarlo simulation method. In the Monte Carlo simulation method, we used thecomputer system R [10].First, we evaluate the probability contents of simplexes. For an integer d ≥ P d and Q d as P d = ( x ∈ R d | x i + √ d ≥ ≤ i ≤ d ) , − x − · · · − x d + √ d ≥ ) ,Q d = ( x ∈ R d | x i − √ d ≥ ≤ i ≤ d ) , − x − · · · − x d + (2 d +1) √ d ≥ ) P d and Q d are simplexes, and they are in general position and bounded.In the Monte Carlo method, we generated 1 , ,
000 simulations from a normaldistribution and computed the fraction of samples that fell into simplexes.The probability contents obtained by the HGM and Monte Carlo methodsare given in Tables 1 and 2. We also show the computational times for the HGMin the tables. d HGM time of HGM(s) MC2 0.285205 0.00 0.28493 0.251995 0.00 0.24934 0.241744 0.01 0.24295 0.242724 0.02 0.24286 0.250219 0.09 0.23947 0.261920 0.32 0.25728 0.276510 1.04 0.27879 0.293138 3.15 0.285910 0.311198 9.51 0.3072Table 1: The probability content of P d as obtained by the HGM and MonteCarlo methods. d HGM time of HGM(s) MC2 5.1758e-02 0.00 5.1917e-023 7.0235e-03 0.00 7.0850e-034 6.3101e-04 0.00 6.0400e-045 3.9722e-05 0.02 5.5000e-056 1.8042e-06 0.10 3.0000e-067 5.9878e-08 0.30 0.0000e+008 1.4799e-09 0.85 0.0000e+009 1.1393e-11 2.25 0.0000e+0010 1.2861e-11 5.74 0.0000e+00Table 2: The probability content of Q d as obtained by the HGM and MonteCarlo methods.Note that the accuracy of the Monte Carlo method is low when the proba-bility content of Q d is very small, and for dimensions greater than 6, the MonteCarlo method could not evaluate the probability content of Q d . The number ofsamples was not enough to evaluate the probability.Next, we estimate the probability content of a simplicial cone. For an integer14 ≥
2, we define a polyhedron C d as C d = ( x ∈ R d | d X i =1 a ij x i + √ d ≥ ≤ j ≤ d ) ) ,a ij = ( i + j ) /
100 ( i < j )1 ( i = j )0 ( i > j ) . We can evaluate the multivariate normal probability of a simplicial cone usingthe method presented in Subsection 5.2. Table 3 shows the probability contentof C d evaluated by the HGM and Monte Carlo methods, and the computa-tional times for the HGM. In the Monte Carlo method, we generated 1 , , References [1] H. Edelsbrunner. The union of balls and its dual shape.
Discrete & Com-putational Geometry , 13:415–440, 1995.[2] H. Hashiguchi, Y. Numata, N. Takayama, and A. Takemura. Holonomicgradient method for the distribution function of the largest root of Wishartmatrix.
Journal of Multivariate Analysis , 117:296–312, 2013.[3] T. Koyama. Holonomic modules associated with multivariate normal prob-abilities of polyhedra. http://arxiv.org/abs/1311.6905 , 2013.[4] T. Koyama and A. Takemura. Calculation of orthant probabilities by theholonomic gradient method. http://arxiv.org/abs/1211.6822 , 2012.155] T. Koyama and A. Takemura. Holonomic gradient method for distribu-tion function of a weighted sum of noncentral chi-square random variables. http://arxiv.org//abs/1503.00378 , 2015.[6] Tamio Koyama, Hiromasa Nakayama, Katsuyoshi Ohara, Tomonari Sei,and Nobuki Takayama. Software packages for holonomic gradient method.In Hoon Hong and Chee Yap, editors,
Mathematical Software – ICMS 2014 ,volume 8592 of
Lecture Notes in Computer Science , pages 706–712. SpringerBerlin Heidelberg, 2014.[7] D. Q. Naiman and H. P. Wynn. Inclusion–exclusion–bonferroni identitiesand inequalities for discrete tube–like problems via euler characteristics.
The Annals of Statistics , 20(1):43–76, 1992.[8] H. Nakayama, K. Nishiyama, M. Noro, K. Ohara, T. Sei, N. Takayama,and A. Takemura. Holonomic gradient descent and its application to theFisher-Bingham integral.
Advances in Applied Mathematics , 47:639–658,2011.[9] T. Oaku. Computation of the characteristic variety and the singular locusof a system of differential equations with polynomial coefficients.
JapanJournal of Industrial and Applied Mathematics , 11:485–497, 1994.[10] R Core Team.
R: A Language and Environment for Statistical Computing .R Foundation for Statistical Computing, Vienna, Austria, 2012. ISBN 3-900051-07-0.[11] Grady Weyenberg, Daniel K Howe, and Ruriko Yoshisa. Nor-malizing kernels in the billera-holmes-vogtmann treespace. http://arxiv.org/pdf/1506.00142v1 , 2015.[12] G. M. Ziegler.