Cumulative Function N Dimensional Gaussians 12.2013

N-DIMENSIONAL CUMULATIVE FUNCTION, AND OTHER USEFUL
FACTS ABOUT GAUSSIANS AND NORMAL DENSITIES
l Bensimhoun (reviewed by David Arnon)

by Michae
Jerusalem, 06/2009
A. 1-d Gaussians and normal densities

The 1-d normal density with mean and standard deviation > 0 is given
by
2
2
1
e(x) /2 .
N (x, , ) =
2
Its variance can be shown to be 2 .

The cumulative function of N (x, 0, 1) cannot be computed explicitly, and is
usually denoted by :
(t) =
N (x, 0, 1)dx.
In the general case, performing the change of variable x = (x )/ in

t
the integral F (t) = N (x, , )dx, the cumulative function of N (x, , ) is
easily seen to be
F (x) = ((x )/).
Notice that F is invertible, with inverse equal to
F 1 (y) = 1 (y) + ,
as can be easily checked.
The Mahalanobis distance of a point x, to the Gaussian that composes the
normal density, is by denition
r=
|x |
;
it measures how many standard deviations is the point x far from the Gaussian.
Given R > 0, let us compute the probability that a point falls at a Mahalanobis distance r R from the Gaussian above. This probability is equal
to
+R
N (x, , )dx,
R
and since a Gaussian is symmetric about its mean, it is equal to

+R
[
]
2
N (x, , )dx = 2 F ( + R) F () = 2F ( + R) 1.
Conversely, if the probability p that a point falls within a certain Mahalanobis distance R from the Gaussian is known, one can compute R explicitly:
[
(p + 1)
]
p+1
p = 2F ( + R) 1 = F ( + R) =
= R = F 1
/
2
2
(p + 1)
= R = 1
.
2
B. 2-d Gaussians and bivariate normal densities
The bivariate normal density with mean = (1 , 2 ) and covariance matrix
is (with z = (x, y))
p(x, y) =
1
1
exp[ (z )1 (z )t ],
2
2 ||
where is assumed to be symmetric and positive denite. It is not dicult to

show that the most general form of such two-dimensional a matrix is
]
[
1 2
12
,
=
1 2
22
where 1 and 2 are dierent from 0, and || < 1. Its inverse matrix 1 can
be shown to be
]
[
1
/(1 2 )
1/12
1
.
=
1/22
1 2 /(1 2 )
The Mahalanobis distance of a point z = (x, y), to the Gaussian that composes the bivariate law above, is by denition given by r such that
r2 = (z )1 (z )t .
Explicitly,
r2 =
[(
( x )( y ) ( y )2 ]
1
x 1 )2
1
2
2
2
+
.
1 2
1
1
2
2
Since || < 1, the locus of the points for which r is constant are ellipses. It
is not dicult to see that a parametrization of such an ellipse at level r2 is
x = r1 cos + 1 ,
y = r2 ( cos + 1 2 sin ) + 2 .
2
Interestingly, if M is the Choleski decomposition of (M is an upper triangular matrix such that M t M = ), then this parametrization can be written
(x, y) = r(cos , sin )M + (1 , 2 );
this can be shown with ease by the explicit computation of M , namely
[
]
1
M=
,
0 2 1 2
or more even simply using the fact that the parametrization above is nothing
else than the locus of the points ruM + , where u = (x, y) is the set of points
such that uut = 1. To see this, it suces to put z = ruM + in the expression
r2 = (z )1 (z )t . This Choleski form of the parametrization holds also
in higher dimensions, and can be obtained eciently in computers.
Another important fact is that the directions of the two axis of the ellipse
above are given by the eigenvectors
of
, while the extent of the great and
small radius are respectively r 1 and r 2 (1 and 2 being respectively the

largest and smallest eigenvalues of ). In place of proving this directly, this
can be proved using the parametrization above. To this end, one can suppose
without loss of generality that (1 , 2 ) = 0. Denoting u = (cos , sin ), we have
seen that the ellipse is the locus of the points z = (x, y) such that
z = ruM,
with
uut = 1.
The great radius of the ellipse is the segment OP , where O is the center of the
ellipse (that is, (0, 0)), and P = (x, y) is one of the two points of the ellipse that
are located at maximal distance from the center of the ellipse. Similarly, the
small radius of the ellipse is given by one of the two points of the ellipse that
are located at minimal distance from the center. In other word, we have to nd
z and u such that uut = 1 and
zz t = x2 + y 2 = r2 uM M t ut = max,
and similarly zz t = r2 uM M t ut = min .
But this is the classic Rayleigh problem, and it is well known that in the maximum case, zz t /r2 = 1 is the largest eigenvalue of M M t , and ut = ut1 is its
corresponding normalized eigenvector. In the minimum case, zz t /r2 = 2 is
the smallest eigenvalue of M M t and ut = ut2 is its corresponding normalized
2
eigenvector.
This gives the desired vectors OP = r u1,2 M , and the extents

r 1 and r 2 of the great and small radii of the ellipse respectively. This
does not prove, yet, that 1,2 are eigenvalues of , since = M t M = M M t .
t = r2 M t ut are in fact eigenvectors of , and
Nevertheless, the vectors OP
1,2
1,2 are their corresponding eigenvalues in : indeed, there holds
M t ut1,2 = M t M M t ut1,2 = M t (1,2 ut1,2 ) = 1,2 M t ut1,2 ;
This shows that 1,2 are the eigenvalues of , with corresponding eigenvectors
proportional to the great and small radii of the ellipse, as was to be shown.
3
Now let p(x, y) be the general bivariate normal density, with as above.
We hope to compute the cumulative function of this Gaussian, as a function
of r. More precisely, our aim is to compute the probability that a point falls
inside the ellipse given by the parametrization above, for a given Mahalanobis
distance r = R.
To this end, let us rst compute the Jacobian of the parametrization above.
It is equal to

1
cos
1 r
sin
x/r x/

=

y/r y/ 2 ( cos + 1 2 sin ) 2 r( sin + 1 2 cos )

= 1 2 1 2 r = r,
as can be easily checked.
To compute the integral of the Gaussian, we can suppose without loss of
generality that = 0 (since integrals are invariant under translation). The
integral to compute is
1
2 ||
dr
r 2
||re 2 .
It is equal to
re
0
r 2
2
dr =
r 2
2
R2 /2
dr /2 =
er dr
= 1 eR
/2
Thus, the desired cumulative function is given by

F (r) = 1 er
/2
It is invertible, an its inverse can be computed explicitly; in other words, if

the probability p that a point falls at a certain Mahalanobis distance r from the
Gaussian is known, then r can be computed explicitly: Put
p = 1 er
/2
Then
r2 /2 = ln(1 p),
or
r=
2 ln(1 p).
Hence,
r = F 1 (p) =
2 ln(1 p).
C. n-dimensional Gaussian and multivariate normal densities

Now, let us examine the n-dimensional case. The n-variate normal density with mean = (1 , 2 , . . . , n ) and covariance matrix is, with z =
(x1 , x2 , . . . , xn ),
p(z) =
1
(2)n/2
||
[
]
exp 12 (z )1 (z )t ,
where is assumed symmetric and positive denite. The Mahalanobis distance

of a point z = (x1 , x2 , . . . , xn ), to the Gaussian that composes the multivariate
normal law above, is by denition given by r such that
r2 = (z )1 (z )t .
The arguments used above in the 2-dimensional case generalize immediately
and show that:
a. The locus of the points z = (x1 , x2 , . . . , xn ) whose Mahalanobis distance
is constant and equal to r is given by
z = ruM + ,
where = M T M is the Choleski decomposition of , and u is the set of points
such that uut = 1 (that is, the unit sphere). This locus can be shown to be a
n-dimensional ellipsoid.
b. The directions of the axis of the ellipsoid above are given by the eigenvectors of , while their respective extents are given by
r
i ,
1 , 2 , . . . , n being the eigenvalues of .

It is more dicult (but not impossible) to generalize the 2-d argument above
in order nd the cumulative function of the multivariate normal law. Instead, we
shall nd another way. Here is a well known fact about normal laws: Let X =
(X1 , X2 , . . . , Xn ) be a n-dimensional random normal variable, with covariance
matrix . Assume furthermore that X is centered (that is, its mean is 0). Then
= E(X T X).
Indeed, if we denote Z = (x1 , x2 , . . . , xn ) and C =
1
,
(2)n/2 ||
we have to
compute the integral

E(X T X) = C
1
Z T Z exp( Z1 Z T )dx1 dx2 dxn .
2
Let = M T M be the Choleski decomposition of . Notice that M is invertible

since
det() = det(M T ) det(M ) = det(M )2 ,
and
det(M ) =
since is positive by hypothesis (i.e. det() > 0). Furthermore,
det(). Let us put Z = Y M , or what is the same, Y = ZM 1 (with Y =
(y1 , . . . , yn )). Thus,
Z1 Z T = Y M M 1 (M T )1 M T Y T = Y Y T .
The Jacobian of the transformation Z = Y M is constant and equal to |M |,
hence the integral above becomes

1
C
M T Y T Y M exp( Y Y T )|M |dy1 dy2 dyn .
2
T
Extracting the
matrices M and M from the integral, and taking into account
that |M | = ; the expression above becomes
(
)

1
1
T
T
T
M
Y Y exp( Y Y )dy1 dy2 dyn M.
2
(2)n/2
The expression inside the parentheses can be shown to be the identity matrix
(sketch of the proof: non-diagonal components cancel by symmetry, and diagonal elements can be computed using the volume of the k-dimensional sphere).
Hence the whole expression is equal to M T IM = M T M = , as was claimed.
Notice that the argument above shows that if a n-dimensional centered
stochastic variable X is normally distributed, with covariance matrix , and
if = M T M is the Choleski decomposition of M , then Y = XM 1 is normally
distributed with covariance matrix equal to I (that is, its gaussian is the product of n independent Gaussians with standard deviation = 1). By the way,
this gives a good mean to implement stochastic generators of the variable X: it
suces to generate n samples y1 , y2 , . . . , yn following separately the 1-d normal
law with standard deviation = 1, and to multiply Y = (y1 , y2 , . . . , yn ) by
M from the right; this gives a n-dimensional stochastic sample of the n-variate
normal law with covariance matrix .
Now, our aim is to nd the cumulative function of the Gaussian which composes X as a function of the Mahalanobis distance r, that is, to nd for each
r, the proportion of instances of X that falls statistically at a Mahalanobis
distance r r from the Gaussian composing X. Without loss of generality,
the variable X can be assumed to be centered. Recall that the Mahalanobis
distance of a point p = (p1 , . . . , pn ) to this Gaussian is given by r2 = p1 pT .
Letting p be an instance of X, we can see r2 itself as as random variable, that is,
r2 = X1 X T , giving rise the variable
r = X1 X T .
Thus, the problem amounts to nd the cumulative function of r. It can be

obtained from the cumulative function of r2 : indeed, if r2 is the cumulative
function of r2 , that is, r2 (c) = P (r2 c), then the cumulative function r of
r is
r (c) = P (r2 c2 ) = r2 (c2 ).
We need a last simplication stage: As above, put X = Y M , where = M T M
is the Choleski decomposition of . We have seen that Y is normal, with
covariance matrix equal to I. Furthermore,
r2 = X1 X T = Y M 1 M T Y T = Y M M 1 (M T )1 M T Y T = Y Y T .
Henceforth, we have only to nd the cumulative function of Y Y T, or in other
words, the cumulative function of the sum of the squares of n independently
normally distributed Gaussian variables, each with mean 0 and standard deviation = 1. Such a distribution is known in the literature as the Chi-square
distribution with n degrees of freedom, and its density has been obtained in
an analytic form: If N1 , . . . , Nn are n independently distributed normal random variables, with mean 0 and standard deviation = 1, then the density
distribution of N12 + + Nn2 is
f (x) =
x(n2)/2 ex/2
,
2n/2 (n/2)
where is the well known Euler Gamma function.

We can now conclude: The cumulative function of the n-variate normal law is
r2 (n2)/2 x/2
x
e
F (r) =
dx.
n/2 (n/2)
2
0
The change of variable x = y 2 , dx = 2y dy, transforms this last expression into
r
2
y n1 ey /2
F (r) =
dy.
(n2)/2 (n/2)
0 2
If n is even, with n = 2m, then since (1) = 1 and (x) = (x 1)(x 1),
1
F (r) = m1
2
(m 1)!
y n1 ey
/2
dy,
and such an integral can be computed explicitly

using integration by parts. If
n is odd, with n = 2m + 1, then since (1/2) = and (x) = (x 1)(x 1),
r
2
1
F (r) =
y n1 ey /2 dy.
m1
2
(m 1/2)(m 3/2) 3/2 2 0
Simplifying by the factors 2 leads to
F (r) =
3 5 7 (n 2) 2
7
r
0
y n1 ey
/2
dy.
By elementary computations, this can be expressed explicitly

by mean of the
x
2
one-dimensional normal cumulative function (x) = 12 ex /2 dx.
Remark: In Matlab, the chi-square density with n-degrees of freedom is given
by the function
chi2pdf(X,n),
while its cumulative distribution function is given by
chi2cdf(s,n).
and its inverse cumulative distribution by
chi2inv(t,n).
Hence, by what has been seen, the cumulative function of any n-variate normal
centered law is
F(r) = chi2cdf(r2 ,n)
(r the Mahalanobis distance to the Gaussian). Conversely, the condence level
of such a law for a given condence c (with 0 c 1) is given by
chi2inv(c, n).
This is the Mahalanobis distance threshold under which a proportion c of points
emitted by any n-variate normal law fall statistically below this threshold.

Cumulative Function N Dimensional Gaussians 12.2013

Uploaded by

Copyright:

Available Formats

Cumulative Function N Dimensional Gaussians 12.2013

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Cumulative Function N Dimensional Gaussians 12.2013

Uploaded by

Copyright:

Available Formats

N-DIMENSIONAL CUMULATIVE FUNCTION, AND OTHER USEFUL

FACTS ABOUT GAUSSIANS AND NORMAL DENSITIES

l Bensimhoun (reviewed by David Arnon)

A. 1-d Gaussians and normal densities

Its variance can be shown to be 2 .

In the general case, performing the change of variable x = (x )/ in

and since a Gaussian is symmetric about its mean, it is equal to

where is assumed to be symmetric and positive denite. It is not dicult to

small radius are respectively r 1 and r 2 (1 and 2 being respectively the

and similarly zz t = r2 uM M t ut = min .

This gives the desired vectors OP = r u1,2 M , and the extents

Thus, the desired cumulative function is given by

It is invertible, an its inverse can be computed explicitly; in other words, if

C. n-dimensional Gaussian and multivariate normal densities

where is assumed symmetric and positive denite. The Mahalanobis distance

1 , 2 , . . . , n being the eigenvalues of .

compute the integral

Let = M T M be the Choleski decomposition of . Notice that M is invertible

Thus, the problem amounts to nd the cumulative function of r. It can be

where is the well known Euler Gamma function.

and such an integral can be computed explicitly

By elementary computations, this can be expressed explicitly

You might also like