Random_Normal Routines
The function random_normal() returns a normally distributed pseudo-random number with zero mean and unit variance.

Unit
QESBPCSRandom

Overloaded Variants
Function Random_Normal: Extended;
Function Random_Normal(RandomGenerator: TRandomGenFunction): Extended;

Declaration
Function Random_Normal: Extended;

Description
Adapted from the following Fortran 77 code ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM. THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.

The algorithm uses the ratio of uniforms method of A.J. Kinderman and J.F. Monahan augmented with quadratic bounding curves.

Parameters
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

function Random_Normal: Extended;
begin
     Result := Random_Normal (DelphiRandom);
End;

Declaration
Function Random_Normal(RandomGenerator: TRandomGenFunction): Extended;

Implementation

function Random_Normal (RandomGenerator: TRandomGenFunction): Extended;
const
     s = 0.449871;
     t = -0.386595;
     a = 0.19600;
     b = 0.25472;
     r1 = 0.27597;
     r2 = 0.27846;
var
     u, v, x, y, q: Extended;
     Done: Boolean;

begin
     //     Generate P = (u,v) uniform in rectangle enclosing acceptance region

     Done := False;
     repeat
          u := RandomGenerator;
          v := RandomGenerator;
          v := 1.7156 * (v - 0.5);

          // Evaluate the quadratic form
          x := u - s;
          Y := abs (v) - t;
          q := Sqr (x) + y * (a * y - b * x);

          // Accept P if inside inner ellipse
          if (q < r1) then
               Done := True
          else if (q <= r2) and (Sqr (v) < -4.0 * Ln (u) * Sqr (u)) then
               Done := True;
     until Done;
     // Return ratio of P's coordinates as the normal deviate
     if u < VSmall then
          raise EMathError.Create (rsDivideByZero);

     Result := v / u;
End;


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