The Mathematica Journal
Departments
Feature Articles
Columns
New Products
New Publications
Classifieds
Calendar
News Bulletins
Mailbox
Letters
Write Us
About the Journal
Staff and Contributors
Submissions
Subscriptions
Advertising
Back Issues
Home
Download this Issue

Bivariate Normal Distribution

Consider computing the cumulative distribution function (CDF) of a multivariate normal (Gaussian) distribution (MultinormalDistribution). Refer to [5, §26.3] or 6. Begin by loading the package Stub.

[Graphics:../Images/index_gr_112.gif]

Next, define a bivariate mean vector

[Graphics:../Images/index_gr_113.gif]

and a ([Graphics:../Images/index_gr_114.gif] symmetric positive definite) covariance matrix.

[Graphics:../Images/index_gr_115.gif]

We can evaluate the CDF at, say,

[Graphics:../Images/index_gr_116.gif]

as follows.

[Graphics:../Images/index_gr_117.gif]
[Graphics:../Images/index_gr_118.gif]

Alternatively, using [5, §26.3.1],

[Graphics:../Images/index_gr_119.gif]

and [5, §26.3.3],

[Graphics:../Images/index_gr_120.gif]

we can compute the numerical integral directly (with reduced PrecisionGoal).

[Graphics:../Images/index_gr_121.gif]
[Graphics:../Images/index_gr_122.gif]
[Graphics:../Images/index_gr_123.gif]

Now consider the following.

[Graphics:../Images/index_gr_124.gif]
[Graphics:../Images/index_gr_125.gif]
[Graphics:../Images/index_gr_126.gif]

This is rather slow, and also small imaginary parts have crept in. Compare this result to that of direct numerical integration.

[Graphics:../Images/index_gr_127.gif]
[Graphics:../Images/index_gr_128.gif]

Using [5, §26.3.20], we can reduce the general bivariate computation to a one-dimensional problem.

[Graphics:../Images/index_gr_129.gif]

All we need to compute is [Graphics:../Images/index_gr_130.gif]. Examples 7--9 of [5, §26.9] demonstrate the computation and use of [Graphics:../Images/index_gr_131.gif]. In particular, example 9 deals with the computation of integrals of the bivariate normal distribution over polygons.

When [Graphics:../Images/index_gr_132.gif], the inner integral can be computed in closed form.

[Graphics:../Images/index_gr_133.gif]

[Graphics:../Images/index_gr_134.gif]
[Graphics:../Images/index_gr_135.gif]
[Graphics:../Images/index_gr_136.gif]

We simplify this expression, noting the domain of [Graphics:../Images/index_gr_137.gif].

[Graphics:../Images/index_gr_138.gif]
[Graphics:../Images/index_gr_139.gif]

Hence [Graphics:../Images/index_gr_140.gif] can be computed as follows.

[Graphics:../Images/index_gr_141.gif]
[Graphics:../Images/index_gr_142.gif]

We need to treat the limiting case [Graphics:../Images/index_gr_143.gif] separately, following [5, §26.3.7-9 and §26.3.15-18].

[Graphics:../Images/index_gr_144.gif]
[Graphics:../Images/index_gr_145.gif]
[Graphics:../Images/index_gr_146.gif]
[Graphics:../Images/index_gr_147.gif]

Note the use of Alternatives for handling both exact ([Graphics:../Images/index_gr_148.gif]) and floating-point values ([Graphics:../Images/index_gr_149.gif]). 

Implementation of §26.3.20 is immediate.

[Graphics:../Images/index_gr_150.gif]

Recomputing the earlier test values is much faster--and we are now working at machine precision.

[Graphics:../Images/index_gr_151.gif]
[Graphics:../Images/index_gr_152.gif]
[Graphics:../Images/index_gr_153.gif]
[Graphics:../Images/index_gr_154.gif]

Another possibility is to convert the integral

[Graphics:../Images/index_gr_155.gif]

to a differential equation,

[Graphics:../Images/index_gr_156.gif]

and then use NDSolve to compute the numerical solution for a range of  [Graphics:../Images/index_gr_157.gif] values for each  [Graphics:../Images/index_gr_158.gif] value. With the initial condition [5, §26.3.19],

[Graphics:../Images/index_gr_159.gif]

we compute the solution for [Graphics:../Images/index_gr_160.gif] using dynamic programming (which can be faster but uses more memory).

[Graphics:../Images/index_gr_161.gif]
[Graphics:../Images/index_gr_162.gif]

We treat the limiting cases [Graphics:../Images/index_gr_163.gif] separately.

[Graphics:../Images/index_gr_164.gif]
[Graphics:../Images/index_gr_165.gif]
[Graphics:../Images/index_gr_166.gif]
[Graphics:../Images/index_gr_167.gif]

To visualize the probability contours, we choose an appropriate PlotRange so as to compare this result with, for example, [5, Figure 26.2].

[Graphics:../Images/index_gr_168.gif]

[Graphics:../Images/index_gr_169.gif]

The agreement is excellent. Here is a surface plot of the cumulative probability.

[Graphics:../Images/index_gr_170.gif]

[Graphics:../Images/index_gr_171.gif]

Returning to the multivariate problem,

[Graphics:../Images/index_gr_172.gif]
[Graphics:../Images/index_gr_173.gif]

we recompute the earlier test values.

[Graphics:../Images/index_gr_174.gif]
[Graphics:../Images/index_gr_175.gif]
[Graphics:../Images/index_gr_176.gif]
[Graphics:../Images/index_gr_177.gif]

We observe another factor of 3 speedup. Note, however, that if any of [Graphics:../Images/index_gr_178.gif] varies then [Graphics:../Images/index_gr_179.gif] changes, so our dynamic code will record a new definition which is not always optimal.


Converted by Mathematica      May 8, 2000

[Tricks of the Trade Index] [Prev Page][Next Page]