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; |
|