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