Stochastic Integration Rules for Infinite Regions.
Stochastic integration rules are derived for infinite integration intervals, generalizing rules developed by Siegel and O'Brien [ SIAM J. Sci. Statist. Comput., 6 (1985), pp. 169--181] for finite intervals. Then random orthogonal transformations of rules for integrals over the surface of the unit m-sphere are used to produce stochastic rules for these integrals. The two types of rules are combined to produce stochastic rules for multidimensional integrals over infinite regions with Normal or Student-t weights. Example results are presented to illustrate the effectiveness of the new rules.
Introduction
A common problem in applied science and statistics is to numerically compute integrals in the form
with . For statistics applications the function p(') may be an unnormalized unimodal
posterior density function and g(') is some function for which an approximate expected value is needed.
We are interested in problems where p(') is approximately multivariate normal (' - Nm (-; \Sigma) or
multivariate Student-t (' - t m (-; \Sigma)). In these cases, a standardizing transformation in the form
can be determined (possibly using numerical optimization), where - is the point where
log(p(')) is maximized, \Sigma is the inverse of the negative of the Hessian matrix for log(p(')) at -, and C
is the lower triangular Cholesky factor for \Sigma The transformed integrals then take the form
w(jjxjj)f(x)dx;
or
Student-t). If the approximation to p(') is good, then f(x)
Accepted for publication in SIAM Journal on Scientific Computing.
y Partially supported by NSF grant DMS-9211640.
can be accurately approximated by a low degree polynomial in x, and this motivates our construction
of stochastic multidimensional polynomial integrating rules for integrals I(f).
This type of integration problem has traditionally been handled using Monte-Carlo algorithms (see
the book by Davis and Rabinowitz, 1984, and the more recent paper by Evans and Swartz, 1992). A
simple Monte-Carlo algorithm for estimating I(f) might use
with the points x i randomly chosen with probability density proportional to w(jjxjj). This Monte-Carlo
algorithm, which is an importance sampling algorithm for the original problem of estimating E(g), is
often effective, but in cases where the resulting f(x) is not approximately constant, the algorithm
can have low accuracy and slow convergence. However, an important feature of simple Monte-Carlo
algorithms is the availability of practical and robust error estimates. If we let oe E denote the standard
error for the sample, then
and
ff
R
\Gammaff
dt.
The new methods that we will describe can be considered a refinement of this Monte-Carlo with importance
sampling algorithm. Simple Monte-Carlo with importance sampling results are exact whenever
the importance modified integrand is constant, but our methods will be exact whenever the importance
modified integrand is a low degree polynomial. Our methods will also provide a robust error estimate
from the sample standard error. The new one-dimensional integration rules that we develop are generalizations
of the rules derived for the interval [-1,1], with weight by Siegel and O'Brien
(1985). Their work extends earlier work by Hammersley and Handscomb (1964), who also considered
the construction of stochastic integration rules for finite intervals. Our work is also partly based on
work by Haber (1969), who introduced the word "stochastic" for generalized Monte-Carlo rules.
Our development of stochastic multidimensional integration rules requires an additional change of
variables to a radial-spherical coordinate system. We let
Z
z t z=1
Z 1w(r)r
z t z=1
w(r)jrj
The numerical approximations to I(f) that we propose to use will be products of stochastic integration
rules for the radial interval (\Gamma1; 1) with weight w(r)jrj m\Gamma1 , and stochastic rules ( of the same polynomial
degree ) for the surface of the unit m-sphere. Averages of properly chosen samples of these rules
will provide unbiased estimates for I(f), and standard errors for the samples can be used to provide
robust error estimates for the I(f) estimates. Our development was partially motivated by the work
of De'ak (1990), who used a transformation to a spherical coordinate system combined with random
orthogonal transformations to develop a method for computing multivariate normal probabilities, but
he did not consider using higher degree rules.
The basic radial integration rules that we use are combinations of the symmetric sums
h(\Gammaae))=2. A radial rule R(h) takes the form
Given points fae i g, the weights fw i g will be determined so that R has polynomial degree 2n + 1. The
points fae i g will be randomly chosen so that R is an unbiased estimate for T
For fixed points fae i g, the selection of the weights is a standard integration rule construction problem.
If we want a degree d rule, it is sufficient that
When k is an odd integer, the equation is automatically satisfied because both R and the integration
operator are symmetric. Define P (h; r) by
Y
Now P (h; r) is a even degree Lagrange interpolating polynomial for h, so it follows from standard
interpolation theory, that P (h;
r)), and the weights fw i g that we need to make R degree are just integrals of the even
degree Lagrange basis functions. We have the following theorem:
Theorem 1 If the points fae i g are distinct non-negative real numbers and the weights fw i g are defined
by
Y
is a degree 2n
We now describe how to choose the points fae i g so that R is an unbiased estimate for T (h). In order
to accomplish this, we need to find a joint probability density function p(ae
Z 1Z 1:::
for any integrable h. We will explicitly show how to do this when and conjecture the
general form for p for n - 3. We will let T use the fact thatZ 1r
The case randomly with density
2ae . Then we have
Z 1C(ae)2ae
Z 1ae
For degree three rule for T (h) is
If we choose ae - 0 randomly with density 2ae m+1 w(ae)
Z 1R 3 (ae) 2ae m+1 w(ae)
dae
Z 1ae
Z 1ae
Z 1ae m+1 w(ae)
Z 1ae
A degree five rule for T (h) is
We will choose ae - 0 and randomly with joint density
where K is determined by the condition
1. We now need to show that EfR 5
(h). There are three terms in R 5 to consider, so we start with the first one, and we find
Z 1Z 1(aeffi)
0:
For the second term we find
Z 1Z 1ae
Z 1Z 1ae
Z 1ae
Z 1ae
Now
Z 1Z 1ae
Z 1Z 1ae
so
Because R 5 is symmetric in ae and ffi , the last term in R 5 also has expected value T (h)=2, so we have
shown that EfR 5 our results in this section with Proposition 1.
Proposition 1 If and the points fae i g, the rules R 2n+1 given by
(1) with weights given by (2), are chosen with probability density proportional to
Y
ae m+1
then R is an unbiased degree 2n
We have proved this for 2. The form for the probability density for n ? 2 is a
conjectured natural generalization of the Siegel and O'Brien Theorem 5.1 (1985). Because of practical
problems associated with generating random ae's from this density when n ? 2 we focus on the
cases.
3 Stochastic Spherical Integration Rules
The spherical surface integrals will be approximated by averages of random rotations of appropriately
chosen rules for the spherical surface. Let
~
with z t
be an integration rule that approximates an integral of a function s(z) over the
surface Um of the unit m-sphere defined by z t z = 1. If Q is an m \Theta m orthogonal matrix then
~
is also an integration rule for s over Um , because Furthermore, if S has polynomial degree
d, then so does SQ , because s(Qz) has the same degree as s(z). If Q is chosen uniformly (see Stewart,
1980) and S has polynomial degree d, then SQ is an unbiased random degree d rule for Um .
There are many choices that could be used for S. We consider rules given in the book by Stroud
(1971, pages 294-296) and the review paper by Mysovskikh (1980, pages 236-237). The rules that we
will combine with radial rules have degree 1, 3 or 5, and we now list them. A simple degree 1 rule is
is the surface content of Um , and z is any point on Um . A simple degree
3 rule is
with the "1" in the j th position. This rule uses 2m values of s(z). A
different degree 3 rule (Mysovskikh, 1980) is
is the j th vertex of a regular m-simplex with vertices on Um . The degree 3 rule -
S 3 is slightly
more expensive to use than S 3 , but it leads to an efficient general degree 5 rule (Mysovskikh, 1980)
(s(\Gammay
The points y j are determined by taking the midpoints of edges of the m-simplex with vertices
projecting those midpoints onto the surface of Um . -
values of s(z). A
degree five rule which extends S 3 (Stroud, 1971, page 294) is
where u j is one of the points in the fully symmetric set that is determined by all possible
permutations and sign changes of the coordinates of the point (r;
2.
4 Stochastic Spherical-Radial Integration Rules
In this section we combine stochastic radial rules with stochastic spherical rules to produce random
rules for I(f). There are many ways that this could be done. A natural approach is to form a stochastic
product rule SR Q;ae (f) from a spherical surface rule S and a radial rule R. Such a rule takes the form
~
If S and R both have degree d, then SR Q;ae (f) will also have degree d (Stroud, 1971, Theorem 2.3-1).
If Q is a uniformly random orthogonal matrix and ae is random chosen with the correct density for R,
then SR Q;ae (f) will be an unbiased estimate for I(f). We have the following theorem:
Theorem 2 If ae is random with density given by Proposition 1, S has degree 2n+1 and Q is an m \Theta m
uniform random orthogonal matrix, then
w(jjxjj)f(x)dx
whenever f is a degree 2n
w(jjxjj)f(x)dx
for any integrable f .
We give three examples of SR rules. A degree one rule constructed from S 1 and R 1 is
ae
Here Q is unnecessary, because uniform random vectors z from Um give unbiased rules. A degree three
rule constructed from S 3 and R 3 is
A degree five rule constructed from -
S 5 and R 5 is
(w
(w
with ~
and ~
ae (f), SR 3
Q;ae and -
respectively. A
sample of one of these rules can be generated, and the sample average used to estimate I(f). The
standard error for the sample can be used to provide an error estimate. For comparison purposes with
the examples in Section 6, we will use SR 0 (f) to denote the one point rule f(z), with the components
of z chosen from Normal(0,1). SR 0 (f) is just the simple Monte-Carlo rule for I(f) with multivariate
normal weight.
5 Implementation Details and Algorithms
In this section we focus on integrals of the form
w(jjxjj)f(x)dx;
. For integrals of this type, we have determined explicit formulas for
the radial rule weights, along with explicit methods for generating the random radial rule points. We
will also discuss the multivariate Student-t weight
We first consider the rule SR 1
ae . In the case
so 2. Therefore
ae
The probability density for ae is proportional to ae degrees of freedom.
It is a standard statistical procedure to generate a random ae with this density (Monahan, 1987). A
standard procedure for generating uniformly random vectors z from Um , consists of first generating
x with components x i random from Normal(0,1) and setting z = x=jjxjj. However, this combined
procedure for generating random vectors aez must be equivalent to just generating random z from
w(jjzjj). Therefore, all we need to do is generate the components z i random
from Normal(0,1), and this is a simpler procedure. We propose the following algorithm for random
degree one rules:
Degree One Spherical-Radial Rule Integration Algorithm
1. Input ffl, m, f and Nmax .
2.
3. Repeat
(a)
(b) Generate a random x with x i - Normal(0,1).
(c) I +D and
4. Output I - I(f), oe
V and N .
The input ffl is an error tolerance, the input Nmax provides a limit on the time for the algorithm, and the
output oe E is the standard error for the integral estimate I. The algorithm computes I and V using a
modified version of a stable one-pass algorithm (Chan and Lewis, 1979). The unscaled sample standard
error oe E will usually be an error bound with approximately 68% certainty. Users of this algorithm
who desire a higher degree of confidence can scale oe E appropriately. For example, a scale factor of 2
increases the certainty level to approximately 95%.
The error estimates obtained by scaling oe E with this algorithm (and the other algorithms in this
section) should be used with caution for low N values. These error estimates are based on the use of
the Central Limit Theorem to infer that the sample averages SR are approximately Normal. A careful
implementation of the algorithms in this section could include an Nmin parameter and/or use a larger
scale factor for oe E for small N values. For large N , a scaled oe E should provide a robust, statistically
sound error estimate, as long as the multivariate normal model adequately represents the tails in the
posterior density. Posterior densities with thicker tails are often more efficiently and reliably handled
using a multivariate Student-t model. One technique for monitoring this is discussed by Monahan and
Genz (1996).
If we consider the Student-t weight, we can see that the density for ae is proportional to r
a change of variable shows this to be proportional to a Beta( m
(see Devroye, 1986, for generating methods), so the random ae's and the uniformly random vectors z
from Um , needed for SR 1 can easily be generated. We can also show jU m jT so the formula for
is the same as the formula for the multivariate Normal case. By making appropriate changes
to line 3(b) and 3(c) of the previous algorithm, a modified algorithm could be produced.
Next, we consider the rule SR 3
Q;ae . Integration by parts with
and therefore
The probability density for ae is proportional to ae m+1 e \Gammaae 2 =2 , a Chi density with m+2 degrees of freedom.
We propose the following algorithm for stochastic degree three rules:
Degree Three Spherical-Radial Rule Integration Algorithm
1. Input ffl, m, f and Nmax .
2.
3. Repeat
(a)
(b) Generate a uniformly random orthogonal m \Theta m matrix Q.
(c) Generate a random ae - Chi(m
(d) For
I +D and
4. Output I - I(f), oe
V and N .
The random orthogonal matrices Q can be generated using a product of appropriately chosen random
reflections (see Stewart, 1980). Other methods are discussed by Devroye (1986, p. 607).
If we consider the Student-t weight case, then integration by parts shows that T
therefore require - ? 2. In this case, SR 3 becomes
Further analysis shows that r m+1
proportional to a Beta( m+2
density, so the random ae's for these SR 3 can easily be generated, and by making appropriate changes
to lines 3(c) and 3(e) of the previous algorithm, a modified algorithm could be produced.
Finally, we consider the rule -
Q;ae . For the weight
and a little algebra shows
In order to develop an algorithm for -
Q;ae , we need a set of regular m-simplex unit vertices g.
We use the set given in Stroud (1971, page 345, correcting a minor misprint), where v
The joint probability density for (ae; ffi ) is proportional to (aeffi) m+1 e \Gamma(ae 2
is not a standard probability density, but there is a transformation to standard densities. Consider the
integral
and make the change of variables
R -=2
Finally, let
sin
sin
Z 1p
The function q m+1 is proportional to a standard Beta(m density. The first
inner integral has the resulting ae ! ffi and the second has ae ? ffi . Because these cases are both
equally likely and -
SR 5 is symmetric in ae and ffi , there is no loss of generality in always using ae ! ffi .
Therefore, we choose r from a Chi(2m and q from a Beta(m density, and then
will be distributed with joint probability density proportional
to (aeffi) m+1 e \Gamma(ae 2 We note here that the same changes of variables could also be
used to provide a practical method for generating the corresponding ae and ffi for the Siegel and O'Brien
(1985) finite interval rules. This question was not addressed in their paper.
We propose the following algorithm for stochastic degree five rules:
Degree Five Spherical-Radial Rule Integration Algorithm
1. Input ffl, m, f and Nmax .
2. compute the m-simplex vertices fv j g.
3. Repeat
(a)
(b) Generate a uniformly random orthogonal m \Theta m matrix Q and set f~v j g.
(c) Generate a random r - Chi(2m
and set
(d) For
ffl For
4. Output I - I(f), oe
V and N .
If we consider the Student-t weight, then it can be shown that T
(- \Gamma2)(- \Gamma4) T 0 , and we must have
4. In this case, we could also produce a formula for -
SR 5 . However, we have not found any easy
method for generating the random ae's and ffi 's needed for R 5 , and so we do not consider this further.
Anyway, for large -, \Gamma( -+m
so the rules that we have already
developed for the multivariate Normal weight should be effective.
A possibly significant overhead cost for the SR 3 and -
SR 5 rules is the generation of the random
orthogonal matrices. Using the algorithm given by Stewart (1980), it can be shown that the cost for
generating one such matrix Q is approximately 4m 3 =3 floating point operations (flops) plus the cost
of generating m 2 =2 Normal(0,1) random numbers. For SR 3 rules the columns of Q are used for the
evaluation points for 2m integrand values, so the overhead cost per integrand value is 2m 2 =3 flops
plus the cost of generating m=4 Normal(0,1) random numbers. Once an integrand evaluation point is
available, we expect the cost for the evaluation of the integrand to be at least O(m), because there
are m components for the input variable for the integrand. However, with application problems in
statistics, the posterior density is often a complicated expression made up of a combination of standard
elementary functions evaluated using the input variable components combined with the problem data
(see the second example in the next section). Therefore, if the O(m) integrand evaluation cost is
measured in flops, we expect the constant in O(m) to be very large, so that the 2m 2 =3 flops for the
generation of the evaluation point for that integrand evaluation should not be significant unless m
is very large. For -
SR 5 rules the Q overhead cost per evaluation point drops to approximately 2m=3
flops (plus the cost of m=4 Normal variates), and this is not significant compared to the integrand
evaluation cost for typical statistics integration problems. We also note here that we need m=2 and
respectively, for the rules SR 1 and SR 0 , per integrand evaluation, so the Normal
variate overhead is higher for the two lowest degree rules. Overall, except for very simple integrands or
large m values, we do not expect the overhead costs for the four rules to be significant compared to the
integrand evaluation cost.
6 Examples
We begin with a simple example, where
The following table of results we obtained using the SR rules:

Table

Test Results from SR Rules
Values I oe E I oe E I oe E I oe E
1000
These results are as expected, with much smaller standard errors for the higher degree rules.
For our second example we use a seven dimensional proportional hazards model problem discussed
by Dellaportas and Wright (1991, 1992). The posterior density is given by
aet
. After we first transform ae using x log(ae), we model p(') with a
multivariate normal approximation. So we use
after computing the mode - and C for log(ae) log(p). We added a scaling constant
e 207:19 , to prevent problems with underflow. In the following table we show results from the use of SR
rules to approximate I(f 2 ), and expected vales for each of the integration variables. The constant S in
the table is a normalizing constant. For each of the respective SR rules we used the computed value of

Table

2: Test Results from SR Rules with 120,000 f 2 Values
Integrand I oe E I oe E I oe E I oe E
For this example, the -
SR 3 and -
SR 5 rule results have standard errors that are smaller than the SR 0
and SR 1 rule standard errors by factors that are on average about one half. Because the decrease in
standard errors is inversely proportional to the square root of the number of samples, approximately
four times as much integrand evaluation work would be needed for this problem when using the SR 0
and SR 1 rules to obtain errors comparable to the errors for the -
SR 3 and -
SR 5 rules. These results are
not as good as those for the previous problem, but the higher degree SR rules are still approximately
four times more efficient than the lower degree rules. The degree five rule was not better than the degree
three rule for this problem. After the standardizing transformation, the problem is apparently close
enough to multivariate normal, so that a rule with degree higher than three does not produce better
results. We did not find any significant difference in running times needed by the four algorithms for
the results in Table 2, and this supports our analysis of the relative importance of overhead costs for
the different rules.
The two examples in this section are meant to illustrate the use of the algorithms given in this
paper. Much more extensive testing is needed in order to carefully compare these algorithms with other
methods available for numerical integration problems in applied statistics. For some of the testing work
that has been recently done with these methods we refer the interested reader to the paper by Monahan
and Genz (1996). Further testing work is still in progress.
7 Concluding Remarks
We have shown how to derive low degree stochastic integration rules for radial integrals with normal
and Student-t weight functions. We have also shown how these new rules can be combined with
stochastic rules for the surface of the sphere, to provide stochastic rules for infinite multivariate regions
with multivariate normal and Student-t weight functions. Results from the examples suggest that
averages of samples of these rules can provide more accurate integral estimates than simpler Monte-Carlo
importance sampling methods. The standard errors from the samples provide robust error estimates
for the new rules.



--R

Computing Standard Deviations: Accuracy
Methods of Numerical Integration

Positive Imbedded Integration in Bayesian Analysis
A Numerical Integration Strategy in Bayesian Analysis

Numerical Quadrature and Cubature
Some Integration Strategies for Problems in Statistical Inference.
Computing Science and Statistics 24
Monte Carlo Methods
Stochastic Quadrature Formulas
Continuous Univariate Distributions-I
An Algorithm for Generating Chi Random Variables

A Comparison of Omnibus Methods for Bayesian Computation
The Approximation of Multiple Integrals by using Interpolatory Cubature Formulae
Unbiased Monte Carlo Integration Methods with Exactness for Low Order Polynomials
The Efficient Generation of Random Orthogonal Matrices with An Application to Condition Estimation
The Approximate Calculation of Multiple Integrals
--TR
multiple integrals;numerical integration;monte carlo;statistical computation