Table Of ContentCalculation of orthant probabilities by the holonomic
gradient method
3
Tamio Koyama∗ and Akimichi Takemura†‡
1
0
2
January 21, 2013
n
a
J
8
Abstract
1
We apply theholonomic gradient method (HGM) introduced by [9] to thecalcu-
]
A lation of orthant probabilities of multivariate normal distribution. The holonomic
N gradient method applied to orthant probabilities is found to be a variant of Plack-
. ett’s recurrence relation ([14]). However an implementation of the method yields
s
c recurrence relations more suitable for numerical computation than Plackett’s re-
[
currence relation. We derive some theoretical results on the holonomic system
2 for the orthant probabilities. These results show that multivariate normal orthant
v
probabilities possess some remarkable properties from the viewpoint of holonomic
2
2 systems. Finally we show that numerical performance of our method is comparable
8
or superior compared to existing methods.
6
.
1
1 1 Introduction
2
1
v: The holonomic gradient method (HGM) introduced by
i Nakayama et al. [9] is a new method of numerical calculation which utilizes algebraic
X
properties of differential equations. This method has many applications in statistics. For
r
a
example, an application to the evaluation of the exact distribution function of the largest
rootofaWishartmatrixisintroducedin[5]andanapplicationtothemaximum likelihood
estimation for the Fisher-Bingham distribution on the d-dimensional sphere is introduced
in [7]. These applications greatly expand the scope of the field of algebraic statistics.
In this paper, we utilize the holonomic gradient method for an accurate evaluation of
the orthant probability
Φ(Σ,µ) = P(X ≥ 0,...,X ≥ 0), (1)
1 d
where the d-dimensional random vector X = (X ,...,X ) ∈ Rd is normally distributed
1 d
with mean µ and covariance matrix Σ, i.e., X ∼ N(µ,Σ). Since evaluation of the orthant
∗Department of Mathematics, Kobe University
†Graduate School of Information Science and Technology, University of Tokyo
‡JST CREST
1
probabilityhasimportantapplicationsinstatisticalpractice, therearemanystudiesabout
it. Genz introduced a method to calculate the orthant probability utilizing the quasi
Monte-Carlo method in [4]. Miwa, Hayter and Kuriki proposed a recursive integration
algorithm to evaluate the orthant probability in [8].
When the mean vector µ is equal to zero, the orthant probability can be interpreted
as the area of a d − 1 dimensional spherical simplex (cf. [1]). In [16], Schla¨fli gave a
classical differential recurrence formula for Φ(Σ,0). Plackett generalized Schla¨fli’s result
and gave a recurrence formula for Φ(Σ,µ) in [14]. He evaluated orthant probabilities
by a recursive integration utilizing his formula. Gassmann implemented the Plackett’s
method in the case of higher dimensions in [3]. Their method is a recursive integration
basedonadifferential recurrence formula, whereas theholonomicgradientmethodutilizes
differential equations. Plackett’s recurrence formula is not suitable for the holonomic
gradient method and in this paper we give a new recurrence formula, which is more
natural from the viewpoint of holonomic gradient method.
Let y = Σ−1µ and x = −1Σ−1. We denote the i-th element of y by y and the (i,j)
2 i
element of x by x . By this transformation of the parameters, the orthant probability in
ij
(1) can be written as
1
(−π)−d/2det(x)1/2exp − ytx−1y g(x,y),
2
(cid:18) (cid:19)
where
∞ ∞ d d
g(x,y) = ... exp x t t + y t dt (dt = dt ...dt ). (2)
ij i j i i 1 d
Z0 Z0 i,j=1 i=1 !
X X
In order to evaluate the orthant probability, it is enough to evaluate g(x,y).
To apply the holonomicgradient method, we need anexplicit formofa Pfaffian system
([7, Section 2]) associated with g(x,y). The Pfaffian system can be obtained from a
holonomic system for g(x,y) (see [9]). For a given d, we can obtain a holonomic system
for g(x,y) by applying the algorithms introduced in [12] with computer algebra systems.
However, for general d, these algorithms can not be applied and we need theoretical
considerations. In this paper, we give a holonomic system for the function g(x,y) and
construct a Pfaffian system from the holonomic system.
Note that the integral g(x,y) satisfies an incomplete A-hypergeometric system ([10])
when the integration domain is not the orthant but a polytope.
The organization of this paper is as follows. In Section 2, we describe a holonomic
system for an integral
d d
g(x,y) = f(t)exp x t t + y t dt, (3)
ij i j i i
Rd !
Z i,j=1 i=1
X X
where f(t) is a holonomic function or a holonomic distribution. In Section 3, we construct
a holonomic system for the function in (2) by utilizing the result in Section 2. Then we
construct a Pfaffian system from the holonomic system. Finally in Section 4 we describe
numerical experiments of the holonomic gradient method.
2
2 Holonomic system associated with the expectation
under multivariate normal distributions
In this section, we consider the holonomic ideal which annihilates the integral (3), which
is the expectation under multivariate normal distributions (for the definition of the holo-
nomic ideal, see [15]). We develop a general theory, where f(t) in (3) is a smooth function
or a distribution in the sense of Schwartz ([17]). This theory is a generalization of the re-
sultintroducedin[6]. InSection3wewillspecializef(t)tobetheindicatorfunctionofthe
positive orthant. Note that the results of this section can be applied to various problems
of the multivariate normal distribution theory, other than the orthant probability.
We denote the ring of differential operators in x with polynomial coefficient by D =
Chx ,...,x ,∂ ,...,∂ i. The operator ∂/∂x is denoted by ∂ . We frequently use the
1 n x1 xn i xi
following rings:
D := Chx ,y ,t ,∂ ,∂ ,∂ : 1 ≤ i ≤ j ≤ d,1 ≤ k ≤ di,
xyt ij k k xij yk tk
D := Chx ,y ,∂ ,∂ : 1 ≤ i ≤ j ≤ d,1 ≤ k ≤ di,
xy ij k xij yk
D := Chx ,∂ : 1 ≤ i ≤ j ≤ di,
x ij xij
D := Cht ,∂ : 1 ≤ i ≤ di. (4)
t i ti
We use the following notation.
d d
1
y = Σ−1µ, x = − Σ−1, h(x,y,t) = x t t + y t . (5)
ij i j i i
2
i,j=1 i=1
X X
2.1 The case of a smooth function
We first consider the case when f(t) is a smooth function. In this case, the integral (3)
converges if the matrix −x is positive definite and f(t) is of exponential growth.
To state the main theorem of this section, we show a general lemma. The notation
f,g,x,y in the following lemma is generic and not related to (5).
Lemma 1. Let D (resp. D ) be the ring of differential operators with polynomial coeffi-
x xy
cient Chx ,∂ : 1 ≤ i ≤ ni (resp. Chx ,y ,∂ ,∂ : 1 ≤ i ≤ n,1 ≤ j ≤ mi). Suppose that
i xi i j xi yj
a holonomic ideal I annihilates a function f(x,y) on Rn×Rm and the function f(x,y) is
rapidly decreasing with respect to the variable y for any x in an open set O ⊂ Rn. Then
the integration ideal of I with respect to the variable y annihilates
g(x) := f(x,y)dy (6)
Rm
Z
defined on the open set O.
Proof. Since the function f(x,y) is rapidly decreasing, the integral of (6) converges. Let
P beinthe integration ideal J = I + m ∂ D ∩D , then we have P•g = P•fdy
j=1 yj x,y x
(cid:16) (cid:17)
P R
3
by the Lebesgue convergence theorem. Since the differential operator P can be written
as
m
P = P + ∂ P ∈ D (P ∈ I, P ∈ D ),
0 yi i x 0 i xy
i=1
X
we have P •fdy = ∂ P •fdy. Since P •f is rapidly decreasing, we have ∂ P •
yi i i yi i
fdy = 0.
R PR R
We now go back to g(x,y) in (3) and use the notation in (4) and (5). Consider a
C-algebra morphism ϕ from D to D defined by
t xy
d
ϕ : D → D t 7→ ∂ , ∂ 7→ −y −2 x ∂ .
t xy i yi ti i ik yk
!
k=1
X
Here, we assume x = x , ∂ = ∂ . Since [ϕ(∂ ),ϕ(∂ )] := ϕ(∂ )ϕ(t )−ϕ(t )ϕ(∂ ) =
ij ji xij xji ti tj ti j j ti
δ , where δ is Kronecker’s delta, ϕ is well-defined as a morphism of C-algebra.
ij ij
Now we state the main theorem in this section.
Theorem 2. Suppose a function f on Rd is smooth and of exponential growth. If dif-
ferential operators P ,...,P ∈ D annihilate f, then the following differential operators
1 s t
annihilate the integral g(x,y) in (3).
ϕ(P ) (1 ≤ k ≤ s), (7)
k
∂ −2∂ ∂ (1 ≤ i < j ≤ d), (8)
xij yi yj
∂ −∂2 (1 ≤ i ≤ d). (9)
xii yi
Moreover, the differential operators (7),(8),(9) generate a holonomic ideal in D if the
xy
differential operators P ,...,P generate a holonomic ideal in D .
1 s t
Proof. For a differential operator
P = c tα1 ···tαd∂β1 ···∂βd ∈ D
α,β 1 d t1 td t
α,β
X
and p ,q ∈ D (i = 1,...,d), we put
i i xyt
P(p ;q ) = c pα1 ···pαdqβ1···qβd. (10)
i i α,β 1 d 1 d
α,β
X
By the assumption, the differential operators
P (ℓ = 1,...,s), ∂ (1 ≤ i ≤ j ≤ d), ∂ (1 ≤ i ≤ d)
ℓ xij yi
annihilate f and generate a holonomic ideal in D . By the lemma introduced in [13,
xyt
Section 3.3], the differential operators
∂h ∂h ∂h
P (t ;∂ − )(ℓ = 1,...,s), ∂ − (1 ≤ i ≤ j ≤ d), ∂ − (1 ≤ i ≤ d)
ℓ i ti ∂t xij ∂x yi ∂y
i ij i
4
annihilate exp(h)f and generate holonomic ideal I in D . As we will show in Lemma
xyt
4 below, the differential operators (7),(8),(9) generate the integration ideal J of I with
respect to the variables t (i = 1,...,d).
i
Since the function exp(h)f is rapidly decreasing, the integration ideal J annihilates
the integral g(x,y) by Lemma 1. Hence, the differential operators (7),(8),
(9) annihilate g(x,y).
Since the integration ideal of a holonomic ideal is also holonomic (see, e.g., [2, Chap
1]), the ideal J is holonomic when the ideal I is holonomic.
Now, we show that the integration ideal J is generated by the differential operators
(7),(8),(9) inthe following two lemmas. Wedenote P ≡ Qwhen P−Q ∈ d D (∂ −
i=1 xyt yi
t ).
i
P
Lemma 3. Let p = ∂ − y − 2 d x t and q = ∂ − y − 2 d x ∂ . For a
i ti i k=i ik i i ti i k=i ik yk
differential operator P ∈ D and a multi index α ∈ Nd = {0,1,2...}d, the following
t P 0 P
equivalence relations hold.
P(t ;p ) ≡ P(∂ ;q ), (11)
i i yi i
tαP(∂ ;q ) ≡ ∂αP(∂ ;q ). (12)
yi i y yi i
Here, we use the notation in (10).
Proof. By the straightforward calculation, we have the following relations for 1 ≤ i,j ≤ d.
[p ,q ] = 0, (13)
i j
[p ,p ] = [q ,q ] = 0, (14)
i j i j
[q ,t ] = [q ,∂ ] = δ . (15)
i j i yj ij
In order to prove (11), it is sufficient to show
tα1 ···tαdpβ1···pβd ≡ ∂α1 ···∂αdqβ1···qβd (α,β ∈ Nd). (16)
1 d 1 d y1 yd 1 d 0
By the induction on β , we have a relation
i
t qβi ≡ ∂ qβi (1 ≤ i,j ≤ d, β ∈ N ). (17)
j i yj i i 0
When β = 0, the relation (17) holds clearly. Suppose that the relation (17) holds for β .
i i
Then we have
t qβi+1 = t q qβi = (q t −δ )qβi
j i j i i i j ij i
≡ (q ∂ −δ )qβi = ∂ q qβi = ∂ qβi+1.
i yj ij i yj i i yj i
By induction, the relation (17) holds for any β .
i
By the relations (13) and (17), we have
pβ1···pβd ≡ qβ1···qβd (β ∈ N ). (18)
1 d 1 d i 0
5
Now we prove the relation (16) by induction on the multi index α ∈ Nd. When α = 0,
0
the relation (16) holds because of (18). Suppose the relation (16) holds for α, then we
have
t tα1 ···tαdpβ1···pβd ≡ t ∂α1 ···∂αdqβ1···qβd
i 1 d 1 d i y1 yd 1 d
= ∂α1 ···∂αd( qβj)t qβi (by (14) and (15))
y1 yd j i i
j6=i
Y
≡ ∂α1 ···∂αd( qβj)∂ qβi (by (17))
y1 yd j yi i
j6=i
Y
≡ ∂ ∂α1 ···∂αdqβ1···qβd (by (14) and (15)).
yi y1 yd 1 d
By induction, the relation (11) holds for any α.
Finally, the relation (12) holds by (14),(15) and (17).
Lemma 4. The integration ideal J is generated by the differential operators
(7),(8),(9).
Proof. Note that the ideal I is generated by differential operators
d
P (t ;∂ −y −2 x t ) (ℓ = 1,...,s),
ℓ i ti i ik i
k=i
X
∂ −2t t (1 ≤ i < j ≤ d),
xij i j
∂ −t2, ∂ −t (1 ≤ i ≤ d).
xii i yi i
By (11), the ideal I is generated by
d
P (∂ ;∂ −y −2 x ∂ ) (ℓ = 1,...,s), (19)
ℓ yi ti i ik yk
k=i
X
∂ −2∂ ∂ (1 ≤ i < j ≤ d), (20)
xij yi yj
∂ −∂2 (1 ≤ i ≤ d), (21)
xii yi
∂ −t (1 ≤ i ≤ d). (22)
yi i
˜ ˜
We denote by J the left ideal generated by (7),(8),(9). Clearly we have J ⊂ J. If P
is a differential operator in J, then P can be written as
d d
P = Q+ ∂ R + S (∂ −t ) (Q ∈ I, R ,S ∈ D ) (23)
ti i i yi i i i xyt
i=1 i=1
X X
by the definition of integration ideal. The differential operator Q is written as a linear
combination of the differential operators (19)–(21) with D coefficients. By the second
xyt
term of the right-hand side of (23), we can assume without loss of generality that the
variables ∂ do not appear in these coefficients. By the relation (12) in Lemma 3, we can
ti
6
assume that the coefficient in Q is an element in D . Since any differential operator in
xy
(19)–(21) is an element of J˜+ d ∂ ·D , we can assume Q ∈ J.
i=1 ti xyt
Consider the equation
P
d d
P −Q− S (∂ −t ) = ∂ R .
i yi i ti i
i=1 i=1
X X
We can assume that the variables ∂ ,...,∂ do not appear in S ,...,S . For example, if
t1 td 1 d
t ∂ is a term of S then we can replace S and R to S −t ∂ −1 and R −t (∂ −t )
1 t1 2 2 1 2 1 t1 1 1 y2 2
respectively. The term t ∂ in S is removed. In the same way, we can remove all terms
1 t1 2
which include ∂ .
ti
Expanding both sides and comparing the coefficients of ∂ , we have
ti
P −Q = S (∂ −t ).
i yi i
X
Theright-handsideofthisequationisanelementoftheleftidealJ′ := d D ·(∂ −t ).
i=1 xyt yi i
Let the weight of t be 1 and that of other variables 0, and consider a term order ≺ with
i
P
this weight. The set {t −∂ |1 ≤ i ≤ d} is a Gr¨obner basis of J′ with the order, so that
i yi
the initial term of P −Q has to divide some t . Since P −Q is in D , we have P −Q = 0.
i xy
Thus P ∈ J.
2.2 The case of a non-smooth function
Next we consider the case when f(t) is not smooth. In this case we consider f(t) as a
distribution in the sense of Schwartz ([17]).
Let Ω be a domain defined by
{(x,y)|−x is positive definite}.
For a tempered distribution f on Rd, we can define a function on Ω as
g(x,y) = hf,exp(h(t,y,x))i. (24)
Sinceexp(h(t,y,x))israpidlydecreasingwithrespecttothevariabletwhen−xispositive
definite, the right-hand side of (24) is finite.
A holonomic system for (24) is given as follows.
Theorem 5. If differential operators P ,...,P ∈ D annihilate a tempered distribution
1 s t
f on Rd, then the differential operators (7),(8),(9) annihilate the function g(x,y) in (24).
Moreover, the differential operators (7),(8),(9) generate a holonomic ideal in D if the
xy
differential operators P ,...,P generate a holonomic ideal in D .
1 s t
Proof. We only need to prove that the differential operators (7),(8),(9) annihilate g(x,y).
Let P = c tα∂β and P∗ = c tβ∂α. We have
ℓ αβ t ℓ αβ t
P Pd
P (∂ ;−y −2 x ∂ )hf,exp(h(t,y,x))i
ℓ yi i ik yk
k=i
X
7
d
= hf,P (∂ ;−y −2 x ∂ )exp(h(t,y,x))i
ℓ yi i ik yk
k=i
X
= hf,P (∂ ;−∂ )exp(h(t,y,x))i
ℓ yi ti
= hf,P∗(−∂ ;∂ )exp(h(t,y,x))i
ℓ ti yi
= hf,P∗(−∂ ;t )exp(h(t,y,x))i
ℓ ti i
= hP (t ;∂ )f,exp(h(t,y,x))i
ℓ i ti
= 0.
3 Holonomic system associated with the orthant prob-
ability
In this section we specialize f(t) of the last section to the indicator function of the positive
orthant and we will construct a Pfaffian system associated with the integral (2) for the
orthant probability.
3.1 Generators of the holonomic ideal
At first, we obtain generators of a holonomic ideal which annihilates (2) by Theorem 5.
Let E be the positive orthant in Rd defined by
t = (t ,...,t ) ∈ Rd|t ≥ 0,(i = 1,...,d) ,
1 d i
(cid:8) (cid:9)
and 1 be the indicator function of E. A holonomic ideal which annihilates 1 is given
E E
as follows.
Lemma 6. The indicator function 1 is annihilated by the following differential operators
E
as a distribution.
t ∂ ,...,t ∂ (25)
1 t1 d td
The differential operators (25) generate a holonomic ideal J in the ring D .
t
Proof. At first, we show that the differential operators t ∂ annihilates the function 1 .
i ti E
It suffices to prove for i = 1. Let ϕ(t) be a rapidly decreasing function on Rd. Then, we
have
∞
− ∂ (t ϕ(t))dt = 0
t1 1 1
Z0
for any t ∈ R, 2 ≤ j ≤ d. Integrating both sides with respect to the variables t ,...,t ,
j 2 d
we have ht ∂ 1 ,ϕi = 0. Therefore, the distribution t ∂ 1 is equal to 0.
1 t1 E 1 t1 E
Next, we show that the left ideal J is holonomic. By the Buchberger’s criterion, we
can show that the set of the differential operators in (25) is a Gr¨obner basis of J with the
8
weight w = (0,1). The characteristic variety (see, e.g. [11], [15]) of J is ch(J) = {(t,ξ) ∈
C2d|t ξ = 0, 1 ≤ i ≤ d}. The variety ch(J) can be decomposed as follows,
i i
{(t,ξ)|t = 0,ξ = 0,i ∈ J,j ∈/ J}.
i j
J⊂{1,...,d}
[
Since the Krull dimension of the each component is d, we have dim(ch(J)) = d and the
ideal J is holonomic.
The function g(x,y)in (2) can bewritten asg(x,y) = h1 (t),exp(h(x,y,t))i.This is a
E
case of Theorem 5 in which the distribution f is 1 . A holonomic ideal which annihilates
E
1 is given in Lemma 6. Hence we have the following theorem for g(x,y) in (2).
E
Theorem 7. Differential operators
d
2 x ∂ ∂ +y ∂ +1 (i = 1,...,d, x = x ), (26)
ik yi yk i yi ij ji
k=1
X
∂ −2∂ ∂ (1 ≤ i < j ≤ d), (27)
xij yi yj
∂ −∂2 (1 ≤ i ≤ d) (28)
xii yi
annihilate the function g(x,y) in (2), and generate a holonomic ideal I in the ring D .
xy
3.2 Differential recurrence formula
Next, we give a Pfaffian system associated with (2). In order to write the Pfaffian system,
we define new notation. For J ⊂ [d] = {1,...,d}, we put
h (x,y,x) = x t t + y t ,
J ij i j k k
i∈J j∈J k∈J
XX X
∞ ∞
g (x,y) = ... exp(h (x,y,t))dt ,
J J J
Z0 Z0
where dt = dt . When theset J is empty, we set g = 1. Forexample, thefunctions
J j∈J j ∅
are written as follows when d = 2.
Q
∞
g (x,y) = exp t2x +2t t x +t2x +y t +y t dt dt
{1,2} 1 11 1 2 12 2 22 1 1 2 2 1 2
Z0
= g(x,y), (cid:0) (cid:1)
∞
g (x,y) = exp t2x +y t dt ,
{1} 1 11 1 1 1
Z0
∞ (cid:0) (cid:1)
g (x,y) = exp t2x +y t dt ,
{2} 2 22 2 2 2
Z0
g = 1. (cid:0) (cid:1)
∅
9
Let J be a subset of [d] and put
d
Q = − y +2 x ∂ (j = 1,...,d), (29)
j j jk yk
!
k=1
X
Q = − y +2 x ∂ (j ∈ J). (30)
j,J j jk yk
!
k∈J
X
Lemma 8. The following equation holds.
Q g = Q g = g (j ∈ J). (31)
j J j,J J J\{j}
Proof. Since g is constant with respect to y (ℓ ∈/ J), we have Q g = Q g . We assume
J ℓ j J j,J J
that J = [d] without loss of generality. Applying Q to the integrand of g , we have
j [d]
d
− y +2 x ∂ exp(h(x,y,t))
j jk yk
!
k=1
X
d
= − y +2 x t exp(h(x,y,t))
j jk k
!
k=1
X
= −∂ exp(h(x,y,t)).
tj
Integrating both sides of the equation from 0 to ∞ with respect to t , we have
j
∞
Q exp(h(x,y,t))dt = exp h (x,y,t) .
j j [d]\{j}
Z0
(cid:0) (cid:1)
Integrating both sides of the equation with respect to the remaining variables, we have
the equation (31).
Remark. By Lemma 8, we have Q Q g = g for i 6= j When the vector y is
i j [d]\{i,j}
equal to zero, this equation can be written as
d
2x +4 x x ∂ +2 x x ∂ g = g .
ij ik jk xkk ik jℓ xkℓ [d]\{i,j}
!
k=1 1≤k6=ℓ≤d
X X
By the transformation of parameters with Σ = (σ ) = −1x−1, we have
kℓ 2
∂
|σ|−1/2g = |σ|−1/2g .
∂σJ [d]\{i,j}
kℓ
(cid:0) (cid:1)
It can be easily checked that this equation corresponds to the classical Schla¨fli’s formula.
For J = {j ,...,j } ⊂ [d], (j < j < ··· < j ), we denote
1 s 1 2 s
x = (x ) , y = (y ,...,y )T,
J jkjℓ 1≤k,ℓ≤s J j1 js
1
Σ = − x−1 = (σJ), µ = Σ y = (µJ ,...,µJ ).
J 2 J ij J J J j1 js
The following differential recurrence formula holds.
10
