| Random_Binomial1 Routines | 
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 | 
| n | Number of Bernoulli Trials, must be positive. | 
| p | 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 FloatsImplementation
 
| 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; | 
|  |