| Random_Binomial2 Routines | 
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 | 
| n | 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 FloatsImplementation
 
| 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; | 
|  |