Random_MVNorm Routines |
Unit
QESBPCSRandom
Overloaded Variants |
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 | H (J) = J'th element of Vector of Means |
d | D (J*(J-1)/2+I) = (I,J)'th element of Variance Matrix (J > = I) |
x | X(J) = J'th element of delivered Vector |
f | 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 FloatsImplementation
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; |
|