| Random_Poisson Routines |
Unit
QESBPCSRandom
| Overloaded Variants |
| Function Random_Poisson(var mu: Extended): Int64; |
| Function Random_Poisson(var mu: Extended; RandomGenerator: TRandomGenFunction): Int64; |
Declaration
Function Random_Poisson(var mu: 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.
For details see: Ahrens, J.H. and Dieter, U. Computer Generation of Poisson Deviates From Modified Normal Distributions. ACM Trans. Math. Software, 8, 2 (June 1982),163-179
| Parameters |
| mu | The mean of the Poisson distribution. |
| 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_Poisson (var mu: Extended): Int64;
begin
Result := Random_Poisson (mu, DelphiRandom);
End; |
Declaration
Function Random_Poisson(var mu: Extended; RandomGenerator: TRandomGenFunction): Int64;
Description
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. GENerate POIsson random deviate Function Generates a single random deviate from a Poisson distribution with mean mu. Arguments mu : The mean of the Poisson distribution from which a random deviate is to be generated. REAL mu Method For details see: Ahrens, J.H. and Dieter, U. Computer Generation of Poisson Deviates From Modified Normal Distributions. ACM Trans. Math. Software, 8, 2 (June 1982),163-179 TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL Implementation
function Random_Poisson (var mu: Extended;
RandomGenerator: TRandomGenFunction): Int64;
const
a0 = -0.5;
a1 = 0.3333333;
a2 = -0.2500068;
a3 = 0.2000118;
a4 = -0.1661269;
a5 = 0.1421878;
a6 = -0.1384794;
a7 = 0.1250060;
var
b1, b2, c, c0, c1, c2, c3, del, difmuk, e, fk, fx, fy, g,
omega, px, py, t, u, v, x, xx: Extended;
s, d, p, q, p0: Extended;
j, k, kflag: Integer;
l, m: Integer;
pp: array [1..35] of Extended;
fact: array [1..10] of Extended;
FirstK0: Boolean;
begin
u := 0;
e := 0;
fk := 0;
difmuk := 0;
fact [1] := 1; // Factorial 0
fact [2] := 1;
fact [3] := 2;
fact [4] := 6;
fact [5] := 24;
fact [6] := 120;
fact [7] := 720;
fact [8] := 5040;
fact [9] := 40320;
fact [10] := 362880; // Factorial 9
Result := 0;
if (mu > 10.0) then
begin
// C A S E A. (RECALCULATION OF S, D, L IF MU HAS CHANGED)
FirstK0 := False;
s := Sqrt (mu);
d := 6.0 * Sqr (mu);
// THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
// PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484)
// IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .
l := Trunc (mu - 1.1484);
// STEP N. NORMAL SAMPLE - random_normal() FOR STANDARD NORMAL DEVIATE
g := mu + s * Random_Normal (RandomGenerator);
if (g > 0.0) then
begin
Result := Trunc (g);
// STEP I. IMMEDIATE ACCEPTANCE IF ival IS LARGE ENOUGH
if (Result >= l) then
Exit;
// STEP S. SQUEEZE ACCEPTANCE - SAMPLE U
fk := Result;
difmuk := mu - fk;
u := RandomGenerator;
if (d * u >= Sqr (difmuk) * difmuk) then
Exit;
end;
// STEP P. PREPARATIONS FOR STEPS Q AND H.
// (RECALCULATIONS OF PARAMETERS IF NECESSARY)
// .3989423=(2*PI)**(-.5) .416667E-1=1./24. .1428571=1./7.
// THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
// APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
// C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.
omega := 0.3989423 / s;
b1 := 0.4166667E-1 / mu;
b2 := 0.3 * Sqr (b1);
c3 := 0.1428571 * b1 * b2;
c2 := b2 - 15.0 * c3;
c1 := b1 - 6.0 * b2 + 45.0 * c3;
c0 := 1.0 - b1 + 3.0 * b2 - 15.0 * c3;
c := 0.1069 / mu;
kflag := -1;
if (g >= 0.0) then
begin
// 'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)
kflag := 0;
FirstK0 := True;
end;
// STEP E. EXPONENTIAL SAMPLE - random_exponential() FOR STANDARD EXPONENTIAL
// DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
// (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)
repeat // 50
if not FirstK0 then
begin
repeat // 50
e := Random_Exponential (RandomGenerator);
u := RandomGenerator;
u := u + u - 1.0;
if u < 0 then
t := 1.8 - abs (e)
else
t := 1.8 + abs (e);
until t > -0.6744;
Result := Trunc (mu + s * t);
fk := Result;
difmuk := mu - fk;
// 'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)
kflag := 1
end
else
FirstK0 := False;
// STEP F. 'SUBROUTINE' F. CALCULATION OF PX, PY, FX, FY.
// CASE ival < 10 USES FACTORIALS FROM TABLE FACT
if (Result < 10) then // 70
begin
px := -mu;
py := XtoY (mu, Result) / fact [Result + 1];
end
else
begin
// CASE ival >= 10 USES POLYNOMIAL APPROXIMATION
// A0-A7 FOR ACCURACY WHEN ADVISABLE
// .8333333E-1=1./12. .3989423=(2*PI)**(-.5)
del := 0.8333333E-1 / fk; // 80
del := del - 4.8 * Sqr (del) * del;
v := difmuk / fk;
if Abs (v) > 0.25 then
px := fk * Ln (1.0 + v) - difmuk - del
else
px := fk * Sqr (v) * (((((((a7 * v + a6) * v + a5) * v + a4)
* v + a3) * v + a2) * v + a1) * v + a0) - del;
py := 0.3989423 / Sqrt (fk);
end;
x := (0.5 - difmuk) / s; // 110
xx := Sqr (x);
fx := -0.5 * xx;
fy := omega * (((c3 * xx + c2) * xx + c1) * xx + c0);
if (kflag <= 0) then
begin
// STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)
if (fy - u * fy <= py * Exp (px - fx)) then // 40
Exit;
end
else
begin
// STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)
if (c * Abs (u) <= py * Exp (px + e) - fy * Exp (fx + e)) then // 60
Exit;
end;
until False;
end
else
begin
// C A S E B. mu < 10
// START NEW TABLE AND CALCULATE P0 IF NECESSARY
m := MaxXY (1, Trunc (mu));
l := 0;
p := Exp (-mu);
q := p;
p0 := p;
// STEP U. UNIFORM SAMPLE FOR INVERSION METHOD
repeat
u := RandomGenerator;
Result := 0;
if (u <= p0) then
Exit;
// STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
// PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
// (0.458=PP(9) FOR MU=10)
if l <> 0 then
begin
j := 1;
if (u > 0.458) then
j := MinXY (l, m);
for k := j to l do
if (u <= pp [k]) then
begin
Result := K;
Exit;
end;
if l = 35 then
Continue;
end;
// STEP C. CREATION OF NEW POISSON PROBABILITIES P
// AND THEIR CUMULATIVES Q=PP(K)
l := l + 1; // 150
for k := l to 35 do
begin
p := p * mu / k;
q := q + p;
pp [k] := q;
if (u <= q) then
begin
Result := K;
Exit;
end;
end;
l := 35;
until False
end;
End; |
|
|