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.

Next, define a bivariate mean vector

and a ( symmetric positive definite) covariance matrix.

We can evaluate the CDF at, say,

as follows.

Alternatively, using [5, §26.3.1],

and [5, §26.3.3],

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

Now consider the following.

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

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

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

When , the inner integral can be computed in closed form.

We simplify this expression, noting the domain of .

Hence  can be computed as follows.

We need to treat the limiting case  separately, following [5, §26.3.7-9 and §26.3.15-18].

Note the use of Alternatives for handling both exact () and floating-point values (). 

Implementation of §26.3.20 is immediate.

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

Another possibility is to convert the integral

to a differential equation,

and then use NDSolve to compute the numerical solution for a range of   values for each   value. With the initial condition [5, §26.3.19],

we compute the solution for  using dynamic programming (which can be faster but uses more memory).

We treat the limiting cases  separately.

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

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

Returning to the multivariate problem,

we recompute the earlier test values.

We observe another factor of 3 speedup. Note, however, that if any of  varies then changes, so our dynamic code will record a new definition which is not always optimal.

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