g05rz sets up a reference vector and generates an array of pseudorandom numbers from a multivariate Normal distribution with mean vector

a

and covariance matrix

C

Syntax

C#
public static void g05rz( int mode, int n, int m, double[] xmu, double[,] c, double[] r, G05..::..G05State g05state, double[,] x, out int ifail )

Visual Basic
Public Shared Sub g05rz ( _ mode As Integer, _ n As Integer, _ m As Integer, _ xmu As Double(), _ c As Double(,), _ r As Double(), _ g05state As G05..::..G05State, _ x As Double(,), _ <OutAttribute> ByRef ifail As Integer _ )

Visual Basic

Public Shared Sub g05rz ( _
	mode As Integer, _
	n As Integer, _
	m As Integer, _
	xmu As Double(), _
	c As Double(,), _
	r As Double(), _
	g05state As G05..::..G05State, _
	x As Double(,), _
	<OutAttribute> ByRef ifail As Integer _
)

Visual C++
public: static void g05rz( int mode, int n, int m, array<double>^ xmu, array<double,2>^ c, array<double>^ r, G05..::..G05State^ g05state, array<double,2>^ x, [OutAttribute] int% ifail )

F#
static member g05rz : mode : int * n : int * m : int * xmu : float[] * c : float[,] * r : float[] * g05state : G05..::..G05State * x : float[,] * ifail : int byref -> unit

Parameters

mode

Type: System..::..Int32

On entry: a code for selecting the operation to be performed by the method.

$mode = 0$: Set up reference vector only.
$mode = 1$: Generate variates using reference vector set up in a prior call to g05rz.
$mode = 2$: Set up reference vector and generate variates.

Constraint:

mode = 0

1

2

n: Type: System..::..Int32
On entry: $n$ , the number of random variates required.

Constraint: $n \geq 0$ .

m: Type: System..::..Int32
On entry: $m$ , the number of dimensions of the distribution.

Constraint: $m > 0$ .

xmu: Type: array<System..::..Double>[]()[][]
An array of size [m]
On entry: $a$ , the vector of means of the distribution.

c: Type: array<System..::..Double,2>[,](,)[,][,]
An array of size [dim1, m]
Note: dim1 must satisfy the constraint: $dim1 \geq m$
On entry: the covariance matrix of the distribution. Only the upper triangle need be set.

Constraint: $C$ must be positive semidefinite to machine precision.

r: Type: array<System..::..Double>[]()[][]
An array of size [lr]
On entry: if $mode = 1$ , the reference vector as set up by g05rz in a previous call with $mode = 0$ or $2$ .
On exit: if $mode = 0$ or $2$ , the reference vector that can be used in subsequent calls to g05rz with $mode = 1$ .

g05state: Type: NagLibrary..::..G05..::..G05State
An Object of type G05.G05State.

x

Type: array<System..::..Double,2>[,](,)[,][,]

An array of size [dim1, _tdx]

Note: dim1 must satisfy the constraint:

dim1 \geq n

if $mode = 1$ or $2$ , $dim1 \geq n$ ;
if $mode = 3$ or $4$ , $dim1 \geq m$ .

On exit: the array of pseudorandom multivariate Normal vectors generated by the method.

Two possible storage orders are available. If

mode = 1

2

then

x [i - 1, j - 1]

holds the

j

th dimension for the

i

th variate. If

mode = 3

4

this ordering is reversed and

x [j - 1, i - 1]

holds the

j

th dimension for the

i

th variate.

ifail: Type: System..::..Int32%
On exit: $ifail = 0$ unless the method detects an error or a warning has been flagged (see [Error Indicators and Warnings]).

Description

When the covariance matrix is nonsingular (i.e., strictly positive definite), the distribution has probability density function

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

where

m

is the number of dimensions,

C

is the covariance matrix,

a

is the vector of means and

x

is the vector of positions.

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.

g05rz decomposes

C

to find such an

L

. It then stores

m

a

and

L

in the reference vector

r

which is used to generate 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 method will work with a singular 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.

One of the initialization methods (G05KFF not in this release) (for a repeatable sequence if computed sequentially) or (G05KGF not in this release) (for a non-repeatable sequence) must be called prior to the first call to g05rz.

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

Error Indicators and Warnings

Errors or warnings detected by the method:

Some error messages may refer to parameters that are dropped from this interface (LDC, LDX) In these cases, an error in another parameter has usually caused an incorrect value to be inferred.

$ifail = 1$: On entry, $mode = 〈value〉$ .
Constraint: $mode = 0$ , $1$ or $2$ .

$ifail = 2$: On entry, $n = 〈value〉$ .
Constraint: $n \geq 0$ .

$ifail = 3$: On entry, $m = 〈value〉$ .
Constraint: $m > 0$ .

$ifail = 5$: On entry, the covariance matrix $C$ is not positive semidefinite to machine precision.

$ifail = 6$: On entry, $ldc = 〈value〉$ and $m = 〈value〉$ .
Constraint: $ldc \geq m$ .

$ifail = 7$: m is not the same as when r was set up in a previous call.
Previous value of $m = 〈value〉$ and $m = 〈value〉$ .

$ifail = 8$: On entry, lr is not large enough, $lr = 〈value〉$ : minimum length required $= 〈value〉$ .

$ifail = 9$: On entry, state vector has been corrupted or not initialized.

$ifail = 11$: On entry, $ldx = 〈value〉$ and $m = 〈value〉$ .
Constraint: $ldx \geq m$ .

On entry, $ldx = 〈value〉$ and $n = 〈value〉$ .
Constraint: $ldx \geq n$ .

$ifail = 12$: On entry, $_tdx = 〈value〉$ and $m = 〈value〉$ .
Constraint: $_sorder \neq 1$ or $_tdx \geq m$ .

$ifail = 22$: On entry, $_tdx = 〈value〉$ and $n = 〈value〉$ .
Constraint: $_sorder = 1$ or $_tdx \geq n$ .

$ifail = -9000$: An error occured, see message report.
$ifail = -6000$: Invalid Parameters $〈value〉$
$ifail = -4000$: Invalid dimension for array $〈value〉$
$ifail = -8000$: Negative dimension for array $〈value〉$
$ifail = -6000$: Invalid Parameters $〈value〉$

Accuracy

Not applicable.

Parallelism and Performance

None.

Further Comments

The time taken by g05rz is of order

n m^{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 covariance matrix that agrees with

C

within machine precision.

Example

Example program (C#): g05rze.cs

Example program data: g05rze.d

Example program results: g05rze.r