Random_Inv_Gauss Routines
Function generates a Random Variate in [0, Infinity) from a reparameterised generalised inverse Gaussian (GIG) Distribution with Density proportional to GIG**(H-1) * EXP(-0.5*B*(GIG+1/GIG)) using a ratio method.

Unit
QESBPCSRandom

Overloaded Variants
Function Random_Inv_Gauss(const h, b: Extended): Extended;
Function Random_Inv_Gauss(const h, b: Extended; RandomGenerator: TRandomGenFunction): Extended;

Declaration
Function Random_Inv_Gauss(const h, b: Extended): Extended;

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

Parameters
Parameter of Distribution, must be positive.
Parameter of Distribution, must be positive.
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_Inv_Gauss (const h, b: Extended): Extended;
begin
     Result := Random_Inv_Gauss (h, b, DelphiRandom);
End;

Declaration
Function Random_Inv_Gauss(const h, b: Extended; RandomGenerator: TRandomGenFunction): Extended;

Implementation

function Random_Inv_Gauss (const h, b: Extended;
     RandomGenerator: TRandomGenFunction): Extended;
var
     ym, xm, r, w, r1, r2, x: Extended;
     a, c, d, e: Extended;
     Done: Boolean;
begin
     if (h < 0.0) or (b <= 0.0) then
          raise EMathError.Create (rsInvalidValue);

     if (h > 0.25 * b * Sqrt (VLarge)) then
          raise EMathError.Create (rsInvalidValue);

     e := Sqr (b);
     d := h + 1.0;
     ym := (-d + Sqrt (Sqr (d) + e)) / b;
     if (ym < VSmall) then
          raise EMathError.Create (rsInvalidValue);

     d := h - 1.0;
     xm := (d + Sqrt (Sqr (d) + e)) / b;
     if (xm < VSmall) then
          raise EMathError.Create (rsInvalidValue);

     d := 0.5 * d;
     e := -0.25 * b;
     r := xm + 1.0 / xm;
     w := xm * ym;
     a := XtoY (w, -0.5 * h) * Sqrt (xm / ym) * Exp (-e * (r - ym - 1.0 / ym));
     if (a < vsmall) then
          raise EMathError.Create (rsInvalidValue);

     c := -d * Ln (xm) - e * r;

     Done := False;
     x := 0.0;
     repeat
          r1 := RandomGenerator;
          if (r1 > 0.0) then
          begin
               r2 := RandomGenerator;
               x := a * r2 / r1;
               if (x > 0.0) then
               begin
                    if (Ln (r1) < d * Ln (x) + e * (x + 1.0 / x) + c) then
                         Done := True
               end;
          end
     until Done;
     Result := x;
End;


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