Random_MVNorm Routines
Routine to generate an N-Variate Random Normal Vector using a Cholesky Decomposition.

Unit
QESBPCSRandom

Overloaded Variants
Procedure Random_MVNorm(const h, d: TESBFloatVector; var f, x: TESBFloatVector; var ier: Integer);
Procedure Random_MVNorm(const h, d: TESBFloatVector; var f, x: TESBFloatVector; var ier: Integer; RandomGenerator: TRandomGenFunction);

Declaration
Procedure Random_MVNorm(const h, d: TESBFloatVector; var f, x: TESBFloatVector; var ier: Integer);

Description
Adapted from Fortran 77 code from the book: Dagpunar, J. 'Principles of random variate generation' Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 N.B. An extra argument, ier, has been added to Dagpunar's routine

Parameters
H (J) = J'th element of Vector of Means
D (J*(J-1)/2+I) = (I,J)'th element of Variance Matrix (J > = I)
X(J) = J'th element of delivered Vector
F((J-1)*(2*N-J)/2+I) = (I,J)'th element of lower triangular decomposition of Variance Matrix (J <= I).
ier 1 if the input covariance matrix is not positive definite, 0 otherwise.
RandomGenerator Optional Function to use for Uniform Random Number Generator. If omitted, Delphi's Random function is used, and if this is done remember to call Randomize if you don't want repeated values.

Category
Arithmetic Routines for Floats

Implementation

procedure Random_MVNorm (const h, d: TESBFloatVector;
     var f, x: TESBFloatVector; var ier: Integer);
begin
     Random_MVNorm (h, d, f, x, ier, DelphiRandom);
End;

Declaration
Procedure Random_MVNorm(const h, d: TESBFloatVector; var f, x: TESBFloatVector; var ier: Integer; RandomGenerator: TRandomGenFunction);

Implementation

procedure Random_MVNorm (const h, d: TESBFloatVector;
     var f, x: TESBFloatVector; var ier: Integer;
     RandomGenerator: TRandomGenFunction);
var
     i, j, m, n, n2, fn: Integer;
     y, v: Extended;
begin
     n := Length (h);
     if (n = 0) then
          raise EMathError.Create (rsVectorEmpty);

     fn := n * (n + 1) div 2;
     if (Length (D) < fn) then
          raise EMathError.Create (rsVectorTooSmall);

     SetLength (f, fn);

     ier := 0;

     n2 := 2 * n;
     if (d [0] <= 0.0) then
     begin
          ier := 1;
          Exit;
     end;

     f [0] := Sqrt (d [0]);
     y := 1.0 / f [0];

     for j := 2 to n do
          f [j - 1] := d [j * (j - 1) div 2] * y;

     for i := 2 to n do
     begin
          v := d [i * (i - 1) div 2 + i - 1];
          for m := 1 to i - 1 do
               v := v - Sqr (f [(m - 1) * (n2 - m) div 2 + i - 1]);

          if (v <= 0.0) then
          begin
               ier := 1;
               Exit;
          end;

          v := Sqrt (v);
          y := 1.0 / v;
          F [(i - 1) * (n2 - i) div 2 + i - 1] := v;

          for j := i + 1 to n do
          begin
               v := d [j * (j - 1) div 2 + i - 1];
               for m := 1 to i - 1 do
               begin
                    v := v - f [(m - 1) * (n2 - m) div 2 + i - 1]
                         * f [(m - 1) * (n2 - m) div 2 + j - 1]
               end;
               f [(i - 1) * (n2 - i) div 2 + j - 1] := v * y;
          end;
     end;

     X := Copy (H, 0, n);
     for j := 1 to n do
     begin
          y := Random_Normal (RandomGenerator);
          for i := j to n do
               x [i - 1] := x [i - 1] + f [(j - 1) * (n2 - j) div 2 + i - 1] * y;
     end;
End;


HTML generated by Time2HELP
http://www.time2help.com