Random_T Routines
Function generates a random variate from a Student T Distribution using Kinderman and Monahan;s Ratio Method.

Unit
QESBPCSRandom

Overloaded Variants
Function Random_T(const DF: Integer): Extended;
Function Random_T(const DF: Integer; RandomGenerator: TRandomGenFunction): Extended;

Declaration
Function Random_T(const DF: Integer): 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
DF Degrees of Freedom - 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_T (const DF: Integer): Extended;
begin
     Result := Random_T (DF, DelphiRandom);
End;

Declaration
Function Random_T(const DF: Integer; RandomGenerator: TRandomGenFunction): Extended;

Implementation

function Random_T (const DF: Integer;
     RandomGenerator: TRandomGenFunction): Extended;
var
     s, c, a, f, g: Extended;
     r, x, v: Extended;
     Done: Boolean;
begin
     if (DF < 1) then
          raise EMathError.Create (rsInvalidValue);

     s := DF;
     c := -0.25 * (s + 1.0);
     a := 4.0 / XtoY (1.0 + 1.0 / s, c);
     f := 16.0 / a;
     if (DF > 1) then
     begin
          g := s - 1.0;
          g := XtoY ((s + 1.0) / g, c) * Sqrt ((s + s) / g);
     end
     else
          g := 1.0;

     Done := False;
     x := 0;
     repeat
          r := RandomGenerator;
          if (r <= 0.0) then
               Continue;
          v := RandomGenerator;
          x := (2.0 * v - 1.0) * g / r;
          v := Sqr (x);
          if (v > 5.0 - a * r) then
          begin
               if (DF >= 3) and (r * (v + 3.0) > f) then
                    Continue;
               if r > XtoY (1.0 + v / s, c) then
                    Continue;
          end;
          Done := True;
     until Done;
     Result := x;
End;


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