NAG Library Function Document

nag_ref_vec_multi_normal (g05eac)

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

+− 9 Example

9.1 Program Text

9.2 Program Data

9.3 Program Results

1 Purpose

nag_ref_vec_multi_normal (g05eac) sets up a reference vector for a multivariate Normal distribution with mean vector

a

and variance-covariance matrix

C

, so that nag_return_multi_normal (g05ezc) may be used to generate pseudorandom vectors.

2 Specification

#include <nag.h>

#include <nagg05.h>

void	nag_ref_vec_multi_normal (const double a[], Integer n, const double c[], Integer tdc, double eps, double *r, NagError fail)

3 Description

When the variance-covariance matrix is non-singular (i.e., strictly positive definite), the distribution has probability density function

f (x) = \sqrt{\frac{|C^{- 1}|}{{(2 π)}^{n}}} \exp \{- {(x - a)}^{T} C^{- 1} (x - a)\}

where

n

is the number of dimensions,

C

is the variance-covariance matrix,

a

is the vector of means and

x

is the vector of positions.

Variance-covariance matrices are symmetric and positive semidefinite. Given such a matrix

C

, there exists a lower triangular matrix

L

such that

{L L}^{T} = C

L

is not unique, if

C

is singular.

nag_ref_vec_multi_normal (g05eac) decomposes

C

to find such an

L

. It then stores

n

a

and

L

in the reference vector

r

for later use by nag_return_multi_normal (g05ezc). nag_return_multi_normal (g05ezc) generates a vector

x

of independent standard Normal pseudorandom numbers. It then returns the vector

a + L x

, which has the required multivariate Normal distribution.

It should be noted that this function will work with a singular variance-covariance matrix

C

, provided

C

is positive semidefinite, despite the fact that the above formula for the probability density function is not valid in that case. Wilkinson (1965) should be consulted if further information is required.

4 References

Knuth D E (1981) The Art of Computer Programming (Volume 2) (2nd Edition) Addison–Wesley

Wilkinson J H (1965) The Algebraic Eigenvalue Problem Oxford University Press, Oxford

5 Arguments

1: a[n] – const doubleInput: On entry: the vector of means, $a$ , of the distribution.
2: n – IntegerInput: On entry: the number of dimensions, $n$ , of the distribution.
Constraint: $n > 0$ .
3: c[ $n \times tdc$ ] – const doubleInput: Note: the $(i, j)$ th element of the matrix $C$ is stored in $c [(i - 1) \times tdc + j - 1]$ .

On entry: the variance-covariance matrix of the distribution. Only the upper triangle need be set.
4: tdc – IntegerInput: On entry: the stride separating matrix column elements in the array c.

Constraint: $tdc \geq n$ .
5: eps – doubleInput: On entry: the maximum error in any element of $C$ , relative to the largest element of $C$ .
Constraint: $0.0 \leq eps \leq 0.1 / n$ .
6: r – double **Output: On exit: reference vector for which memory will be allocated internally. This reference vector will subsequently be used by nag_return_multi_normal (g05ezc). If no memory is allocated to r (e.g., when an input error is detected) then r will be NULL on return, otherwise you should use the NAG macro NAG_FREE to free the storage allocated by r when it is no longer of use.
7: fail – NagError *Input/Output: The NAG error argument (see Section 3.6 in the Essential Introduction).

6 Error Indicators and Warnings

NE_2_INT_ARG_LT: On entry, $tdc = 〈value〉$ while $n = 〈value〉$ . These arguments must satisfy $tdc \geq n$ .
NE_2_REAL_ARG_GT: On entry, $eps = 〈value〉$ while 0.1/ $n = 〈value〉$ . These arguments must satisfy $eps \leq 0.1$ /n.
NE_ALLOC_FAIL: Dynamic memory allocation failed.
NE_INT_ARG_LT: On entry, $n = 〈value〉$ .
Constraint: $n \geq 1$ .
NE_NOT_POS_SEM_DEF: Matrix $C$ is not positive semidefinite.
NE_REAL_ARG_LT: On entry, eps must not be less than 0.0: $eps = 〈value〉$ .

7 Accuracy

The maximum absolute error in

{L L}^{T}

, and hence in the variance-covariance matrix of the resulting vectors, is less than

(n \times \max (eps, ε) + (n + 3) ε / 2)

times the maximum element of

C

, where

ε

is the machine precision. Under normal circumstances, the above will be small compared to sampling error.

8 Further Comments

The time taken by nag_ref_vec_multi_normal (g05eac) is of order

n^{3}

It is recommended that the diagonal elements of

C

should not differ too widely in order of magnitude. This may be achieved by scaling the variables if necessary. The actual matrix decomposed is

C + E = {L L}^{T}

, where

E

is a diagonal matrix with small positive diagonal elements. This ensures that, even when

C

is singular, or nearly singular, the Cholesky factor

L

corresponds to a positive definite variance-covariance matrix that agrees with

C

within a tolerance determined by eps.

9 Example

The example program prints five pseudorandom observations from a bivariate Normal distribution with means vector

[\begin{matrix} 1.0 \\ 2.0 \end{matrix}]

and variance-covariance matrix

[\begin{matrix} 2.0 & 1.0 \\ 1.0 & 3.0 \end{matrix}],

generated by nag_ref_vec_multi_normal (g05eac) and nag_return_multi_normal (g05ezc) after initialization by nag_random_init_repeatable (g05cbc).

9.1 Program Text

Program Text (g05eace.c)

9.2 Program Data

None.

9.3 Program Results

Program Results (g05eace.r)

nag_ref_vec_multi_normal (g05eac) (PDF version)

g05 Chapter Contents

g05 Chapter Introduction

NAG C Library Manual

NAG Library Function Documentnag_ref_vec_multi_normal (g05eac)