Random_Binomial2 Routines
Generates a single random deviate from a binomial distribution whose number of trials is N and whose probability of an event in each trial is P.

Unit
QESBPCSRandom

Overloaded Variants
Function Random_Binomial2(const n: Int64; const pp: Extended): Int64;
Function Random_Binomial2(const n: Int64; const pp: Extended; RandomGenerator: TRandomGenFunction): Int64;

Declaration
Function Random_Binomial2(const n: Int64; const pp: Extended): Int64;

Description
Translated to Fortran 90 by Alan Miller from: RANLIB, Library of Fortran Routines for Random Number Generation. Compiled and Written by:

Barry W. Brown

James Lovato

Department of Biomathematics, Box 237

The University of Texas, M.D. Anderson Cancer Center

1515 Holcombe Boulevard

Houston, TX 77030

This work was supported by grant CA-16672 from the National Cancer Institute.

Parameters
Number of Trials, must be positive.
pp 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_Binomial2 (const n: Int64; const pp: Extended): Int64;
begin
     Result := Random_Binomial2 (n, pp);
End;

Declaration
Function Random_Binomial2(const n: Int64; const pp: Extended; RandomGenerator: TRandomGenFunction): Int64;

Implementation

function Random_Binomial2 (const n: Int64; const pp: Extended;
     RandomGenerator: TRandomGenFunction): Int64;
var
     alv, aMaxXYp, f, f1, f2, u, v, w, w2, x, x1, x2, ynorm, z, z2: Extended;
     i: Integer;
     ix, ix1, k, mp: Int64;
     m: Int64;
     p, q, xnp, ffm, fm, xnpq, p1, xm, xl, xr, c, al, xll,
          xlr, p2, p3, p4, qn, r, g: Extended;
     Done: Boolean;
begin
     // SETUP
     p := MinXY (pp, 1.0 - pp);
     q := 1.0 - p;
     xnp := n * p;

     if (xnp > 30.0) then
     begin
          ffm := xnp + p;
          m := Trunc (ffm);
          fm := m;
          xnpq := xnp * q;
          p1 := INT (2.195 * Sqrt (xnpq) - 4.6 * q) + 0.5;
          xm := fm + 0.5;
          xl := xm - p1;
          xr := xm + p1;
          c := 0.134 + 20.5 / (15.3 + fm);
          al := (ffm - xl) / (ffm - xl * p);
          xll := al * (1.0 + 0.5 * al);
          al := (xr - ffm) / (xr * q);
          xlr := al * (1.0 + 0.5 * al);
          p2 := p1 * (1.0 + c + c);
          p3 := p2 + c / xll;
          p4 := p3 + c / xlr;

          // GENERATE VARIATE, Binomial mean at least 30.

          ix := 0;
          repeat;
               u := RandomGenerator; //20
               u := u * p4;
               v := RandomGenerator;

               // TRIANGULAR REGION
               if (u <= p1) then
               begin
                    ix := Trunc (xm - p1 * v + u);
                    Break;
               end;

               // PARALLELOGRAM REGION
               if (u <= p2) then
               begin
                    x := xl + (u - p1) / c;
                    v := v * c + 1.0 - Abs (xm - x) / p1;
                    if (v > 1.0) or (v <= 0) then
                         Continue;
                    ix := Trunc (x);
               end
               else
               begin
                    // LEFT TAIL
                    if (u <= p3) then
                    begin
                         ix := Trunc (xl + Ln (v) / xll);
                         if (ix < 0) then
                              Continue;
                         v := v * (u - p2) * xll
                    end
                         // RIGHT TAIL
                    else
                    begin
                         ix := Trunc (xr - Ln (v) / xlr);
                         if (ix > n) then
                              Continue;
                         v := v * (u - p3) * xlr;
                    end;
               end;

               // DETERMinXYE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST

               k := Abs (ix - m);
               if (k <= 20) or (k >= xnpq / 2 - 1) then
               begin
                    // EXPLICIT EVALUATION
                    f := 1.0;
                    r := p / q;
                    g := (n + 1) * r;
                    if (m < ix) then
                    begin
                         mp := m + 1;
                         for i := mp to ix do
                              f := f * (g / i - r);
                    end
                    else if (m > ix) then
                    begin
                         ix1 := ix + 1;
                         for i := ix1 to m do
                              f := f / (g / i - r);
                    end;

                    if (v > f) then
                         Continue
                    else
                         Break
               end;

               // SQUEEZING USING UPPER AND LOWER BOUNDS ON LOG(F(X))

               aMaxXYp := (k / xnpq) * ((k * (k / 3.0 + 0.625)
                    + 0.1666666666666) / xnpq + 0.5);
               ynorm := -Sqr (k) / (2.0 * xnpq);
               alv := Ln (v);
               if (alv < ynorm - aMaxXYp) then
                    Break;
               if (alv > ynorm + aMaxXYp) then
                    Continue;
               // STIRLING'S (actually de Moivre's) FORMULA TO MACHINE ACCURACY FOR
               //  THE FINAL ACCEPTANCE/REJECTION TEST

               x1 := ix + 1;
               f1 := fm + 1.0;
               z := n + 1 - fm;
               w := n - ix + 1.0;
               z2 := Sqr (z);
               x2 := Sqr (x1);
               f2 := Sqr (f1);
               w2 := Sqr (w);
               if (alv - (xm * Ln (f1 / x1) + (n - m + 0.5) * Ln (z / w)
                    + (ix - m) * Ln (w * p / (x1 * q))
                    + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / f2) / f2) / f2) / f2) / f1 / 166320.0
                    + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / z2) / z2) / z2) / z2) / z / 166320.0
                    + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / x2) / x2) / x2) / x2) / x1 / 166320.0
                    + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / w2) / w2) / w2) / w2) / w / 166320.0)
                    > 0.0) then
               begin
                    Continue;
               end
               else
                    Break;
          until False;
     end
     else
     begin
          // INVERSE CDF LOGIC FOR MEAN LESS THAN 30
          qn := XtoY (q, n);
          r := p / q;
          g := r * (n + 1);

          Done := False;
          repeat
               ix := 0;
               f := qn;
               u := RandomGenerator;
               while (u >= f) do
               begin
                    if (ix > 110) then
                         Done := False
                    else
                    begin
                         Done := True;
                         u := u - f;
                         ix := ix + 1;
                         f := f * (g / ix - r);
                    end;
               end;
          until Done;
     end;

     if (pp > 0.5) then
          ix := n - ix;

     Result := ix;
End;


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