Random_Beta Routines
Function generates a Random variate in [0,1] from a Beta Distribution with density proportional to BETA**(AA-1) * (1-BETA)**(BB-1).

Unit
QESBPCSRandom

Overloaded Variants
Function Random_Beta(const aa, bb: Extended): Extended;
Function Random_Beta(const aa, bb: Extended; RandomGenerator: TRandomGenFunction): Extended;

Declaration
Function Random_Beta(const aa, bb: Extended): Extended;

Description
Uses Cheng'S log logistic method. Adapted from Fortran 77 code from the book: Dagpunar, J. 'Principles of random variate generation' Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9

Function generates a Random variate in [0,1] from a Beta Distribution with density proportional to BETA**(AA-1) * (1-BETA)**(BB-1). Uses Cheng'S log logistic method.

Parameters
aa Shape Parameter for Beta Distribution - must be positive.
bb Shape Parameter for Beta 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_Beta (const aa, bb: Extended): Extended;
begin
     Result := Random_Beta (aa, bb, DelphiRandom);
End;

Declaration
Function Random_Beta(const aa, bb: Extended; RandomGenerator: TRandomGenFunction): Extended;

Implementation

function Random_Beta (const aa, bb: Extended;
     RandomGenerator: TRandomGenFunction): Extended;
const
     aln4 = 1.3862943611198906;
var
     a, b, g, r, s, x, y, z: Extended;
     d, f, h, t, c: Extended;
     Swap: Boolean;
     Done: Boolean;
begin
     if (aa <= 0.0) or (bb <= 0.0) then
          raise EMathError.Create (rsInvalidValue);

     a := aa;
     b := bb;
     swap := b > a;
     if swap then
     begin
          g := b;
          b := a;
          a := g;
     end;
     d := a / b;
     f := a + b;
     if (b > 1.0) then
     begin
          h := Sqrt ((2.0 * a * b - f) / (f - 2.0));
          t := 1.0;
     end
     else
     begin
          h := b;
          t := 1.0 / (1.0 + XtoY (a / (VLarge * b), b));
     end;
     c := a + h;

     Done := False;
     Result := 0.0;
     repeat
          r := RandomGenerator;
          x := RandomGenerator;
          s := Sqr (r) * x;
          if (r < VSmall) or (s <= 0.0) then
               Continue;

          if (r < t) then
          begin
               x := Ln (r / (1.0 - r)) / h;
               y := d * Exp (x);
               z := c * x + f * Ln ((1.0 + d) / (1.0 + y)) - aln4;
               if (s - 1.0 > z) then
               begin
                    if (s - s * z > 1.0) then
                         Continue;
                    if (Ln (s) > z) then
                         Continue;
               end;
               Result := y / (1.0 + y);
          end
          else
          begin
               if 4.0 * s > XtoY (1.0 + 1.0 / d, f) then
                    Continue;
               Result := 1.0
          end;
          Done := True
     until Done;

     if Swap then
          Result := 1 - Result;
End;


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