g05zs:: Random Number Generators (NAG Toolbox)

A two-dimensional random field

Z (x)

ℝ^{2}

is a function which is random at every point

x \in ℝ^{2}

, so

Z (x)

is a random variable for each

x

. The random field has a mean function

μ (x) = 𝔼 [Z (x)]

and a symmetric positive semidefinite covariance function

C (x, y) = 𝔼 [(Z (x) - μ (x)) (Z (y) - μ (y))]

Z (x)

is a Gaussian random field if for any choice of

n \in ℕ

and

x_{1}, \dots, x_{n} \in ℝ^{2}

, the random vector

{[Z (x_{1}), \dots, Z (x_{n})]}^{T}

follows a multivariate Normal distribution, which would have a mean vector

\tilde{μ}

with entries

{\tilde{μ}}_{i} = μ (x_{i})

and a covariance matrix

\tilde{C}

with entries

{\tilde{C}}_{i j} = C (x_{i}, x_{j})

. A Gaussian random field

Z (x)

is stationary if

μ (x)

is constant for all

x \in ℝ^{2}

and

C (x, y) = C (x + a, y + a)

for all

x, y, a \in ℝ^{2}

and hence we can express the covariance function

C (x, y)

as a function

γ

of one variable:

C (x, y) = γ (x - y)

γ

is known as a variogram (or more correctly, a semivariogram) and includes the multiplicative factor

σ^{2}

representing the variance such that

γ (0) = σ^{2}

The functions nag_rand_field_2d_user_setup (g05zq) or nag_rand_field_2d_predef_setup (g05zr) along with nag_rand_field_2d_generate (g05zs) are used to simulate a two-dimensional stationary Gaussian random field, with mean function zero and variogram

γ (x)

, over a domain

[x_{\min}, x_{\max}] \times [y_{\min}, y_{\max}]

, using an equally spaced set of

N_{1} \times N_{2}

points;

N_{1}

points in the

x

-direction and

N_{2}

points in the

y

-direction. The problem reduces to sampling a Gaussian random vector

X

of size

N_{1} \times N_{2}

, with mean vector zero and a symmetric covariance matrix

A

, which is an

N_{2}

N_{2}

block Toeplitz matrix with Toeplitz blocks of size

N_{1}

N_{1}

. Since

A

is in general expensive to factorize, a technique known as the circulant embedding method is used.

A

is embedded into a larger, symmetric matrix

B

, which is an

M_{2}

M_{2}

block circulant matrix with circulant bocks of size

M_{1}

M_{1}

, where

M_{1} \geq 2 (N_{1} - 1)

and

M_{2} \geq 2 (N_{2} - 1)

B

can now be factorized as

B = W Λ W^{*} = R^{*} R

, where

W

is the two-dimensional Fourier matrix (

W^{*}

is the complex conjugate of

W

Λ

is the diagonal matrix containing the eigenvalues of

B

and

R = Λ^{\frac{1}{2}} W^{*}

B

is known as the embedding matrix. The eigenvalues can be calculated by performing a discrete Fourier transform of the first row (or column) of

B

and multiplying by

M_{1} \times M_{2}

, and so only the first row (or column) of

B

is needed – the whole matrix does not need to be formed.

The symmetry of

A

as a block matrix, and the symmetry of each block of

A

, depends on whether the covariance function

γ

is even or not.

γ

is even if

γ (x) = γ (- x)

for all

x \in ℝ^{2}

, and uneven otherwise (in higher dimensions,

γ

can be even in some coordinates and uneven in others, but in two dimensions

γ

is either even in both coordinates or uneven in both coordinates). If

γ

is even then

A

is a symmetric block matrix and has symmetric blocks; if

γ

is uneven then

A

is not a symmetric block matrix and has non-symmetric blocks. In the uneven case,

M_{1}

and

M_{2}

are set to be odd in order to guarantee symmetry in

B

As long as all of the values of

Λ

are non-negative (i.e.,

B

is positive semidefinite),

B

is a covariance matrix for a random vector

Y

which has

M_{2}

‘blocks’ of size

M_{1}

. Two samples of

Y

can now be simulated from the real and imaginary parts of

R^{*} (U + i V)

, where

U

and

V

have elements from the standard Normal distribution. Since

R^{*} (U + i V) = W Λ^{\frac{1}{2}} (U + i V)

, this calculation can be done using a discrete Fourier transform of the vector

Λ^{\frac{1}{2}} (U + i V)

. Two samples of the random vector

X

can now be recovered by taking the first

N_{1}

elements of the first

N_{2}

blocks of each sample of

Y

– because the original covariance matrix

A

is embedded in

B

X

will have the correct distribution.

References

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

Error Indicators and Warnings

Accuracy

Further Comments

Example

function g05zs_example


fprintf('g05zs example results\n\n');

% Use symmetric stable variogram
icov2 = int64(1); 
params = [0.1; 0.15; 1.2];

% Random Field variance
var = 0.5;
% Domain endpoints
xmin = -1;
xmax =  1;
ymin = -0.5;
ymax =  0.5;
% Number of sample points in x and y
ns = [int64(5), 5];
% maximum dimensions for circulant matrix
maxm = [int64(64), 64];
% Scaling factor, rho = 1.
icorr = int64(2);

% Get square roots of the eigenvalues of the embedding matrix
[lam, xx, yy, m, approx, rho, icount, eig, ifail] = ...
  g05zr( ...
         ns, xmin, xmax, ymin, ymax, maxm, var, ...
         icov2, params, 'icorr', icorr);

fprintf('\nSize of embedding matrix = %d\n\n', m(1)*m(2));

% Display approximation information if approximation used
if approx == 1
  fprintf('Approximation required\n\n');
  fprintf('rho = %10.5f\n', rho);
  fprintf('eig = %10.5f%10.5f%10.5f\n', eig(1:3));
  fprintf('icount = %d\n', icount);
else
  fprintf('Approximation not required\n\n');
end

% Initialize state array
genid = int64(1);
subid = int64(1);
seed  = [int64(14965)];
[state, ifail] = g05kf( ...
                        genid, subid, seed);

% Compute s random field realisations
s = int64(5);
[state, z, ifail] = g05zs( ...
                           ns, s, m, lam, rho, state);

% Display realisations

% Set row labels to grid points (column label is realisation number).
rlabs = cell(ns(1)*ns(2), 1);
for j=1:ns(2)
  for i=1:ns(1)
    if i == 1
      rlabs{(j-1)*ns(1)+i} = sprintf('%6.1f%6.1f', xx(i), yy(j));
    else
      rlabs{(j-1)*ns(1)+i} = sprintf('%6.1f     .', xx(i));
    end
  end
end

% Matrix printing parameters
mtitle = 'Random field realisations (x,y coordinates first):';
matrix = 'General';
diag   = 'Non-unit';
fmt    = 'f10.5';
rlabel = 'Character';
clabel = 'Integer';
clabs  = {' '};
ncols  = int64(80);
indent = int64(0);

[ifail] = x04cb( ...
                 matrix, diag, z, fmt, mtitle, rlabel, rlabs, clabel, ...
                 clabs, ncols, indent);

g05zs example results


Size of embedding matrix = 64

Approximation not required

 Random field realisations (x,y coordinates first):
                     1         2         3         4         5
 -0.8  -0.4   -0.61951  -0.93149  -0.32975  -0.51201   1.38877
 -0.4     .    0.74779   1.33518  -0.51237   0.26595   0.30051
  0.0     .   -0.30579   0.51819   0.50961   0.10379   0.36815
  0.4     .    0.53797  -0.53992  -0.86589  -0.37098   0.21571
  0.8     .   -0.61221  -1.04262   0.00007  -1.22614  -0.06650
 -0.8  -0.2    0.01853   0.64126  -0.42978  -0.79178  -0.55728
 -0.4     .   -0.77912   0.81079  -0.60613   0.07280   1.61511
  0.0     .   -0.23198   1.48744  -0.78145   0.10347   0.07053
  0.4     .    0.32356   0.58676   0.05846   0.34828   1.40522
  0.8     .   -1.24085  -0.92512   0.27247  -0.66965   0.67073
 -0.8   0.0   -1.18183  -0.99775   0.03888   0.01789  -0.65746
 -0.4     .    0.26155  -0.01734  -0.14924   0.28886   0.25940
  0.0     .    1.14960   0.48850  -0.59023   0.22795  -0.60773
  0.4     .   -0.32684  -0.09616  -0.63497  -1.06753  -0.64594
  0.8     .    0.10064   1.06148   0.15020  -0.53168  -0.29251
 -0.8   0.2   -1.30595  -0.03899  -0.35549  -0.20589  -0.35956
 -0.4     .   -0.01776   0.84501   0.20406   0.89039  -0.58338
  0.0     .    0.41898   0.93435  -1.10725   0.76913  -0.74579
  0.4     .   -1.37738   1.72404  -0.20558  -1.41877   1.21816
  0.8     .    0.77866   0.84922  -0.65055   0.83518  -0.26425
 -0.8   0.4   -0.65163   0.50492  -0.52463  -1.12816   1.12817
 -0.4     .    0.15437   0.20739  -0.12675   1.27782  -0.26157
  0.0     .    0.20324   0.54670  -1.73909   0.61580   0.17551
  0.4     .   -1.09470   0.83967   0.70226  -0.34259   0.29368
  0.8     .    1.08452   1.23097  -0.36003   1.06884   0.23594

NAG Toolbox: nag_rand_field_2d_generate (g05zs)

▸▿ Contents

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

Error Indicators and Warnings

Accuracy

Further Comments

Example