Random_Binomial1 Routines
Function generates a Random Binomial Variate using C D Kemp's method.

Unit
QESBPCSRandom

Overloaded Variants
Function Random_Binomial1(const n: Integer; const p: Extended): Int64;
Function Random_Binomial1(const n: Integer; const p: Extended; RandomGenerator: TRandomGenFunction): Int64;

Declaration
Function Random_Binomial1(const n: Integer; const p: Extended): Int64;

Description
This algorithm is suitable when many random variates are required with the SAME parameter values for n & p. Reference: Kemp, C.D. (1986). `A modal method for generating binomial variables', Commun. Statist. - Theor. Meth. 15(3), 805-813.

Parameters
Number of Bernoulli Trials, must be positive.
Bernoulli Probability of Success must be in [0, 1].
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_Binomial1 (const n: Integer; const p: Extended): Int64;
begin
     Result := Random_Binomial1 (n, p, DelphiRandom);
End;

Declaration
Function Random_Binomial1(const n: Integer; const p: Extended; RandomGenerator: TRandomGenFunction): Int64;

Implementation

function Random_Binomial1 (const n: Integer; const p: Extended;
     RandomGenerator: TRandomGenFunction): Int64;
var
     ru, rd: Int64;
     r0: Int64;
     u, pd, pu: Real;
     odds_ratio, p_r: Real;
begin
     r0 := Trunc ((n + 1) * p);
     p_r := bin_prob (n, p, r0);
     if p < 1 then
          odds_ratio := p / (1.0 - p)
     else
          odds_ratio := VLarge;

     u := RandomGenerator;
     u := u - p_r;
     if (u < 0.0) then
     begin
          Result := r0;
          Exit;
     end;

     pu := p_r;
     ru := r0;
     pd := p_r;
     rd := r0;

     repeat
          Dec (rd);
          if (rd >= 0) then
          begin
               pd := pd * (rd + 1.0) / (odds_ratio * (n - rd));
               u := u - pd;
               if (u < 0.0) then
               begin
                    Result := rd;
                    Exit;
               end;
          end;

          Inc (ru);
          if (ru <= n) then
          begin
               pu := pu * (n - ru + 1.0) * odds_ratio / ru;
               u := u - pu;
               if (u < 0.0) then
               begin
                    Result := ru;
                    Exit;
               end;
          end;
     until False;
End;


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