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