Random_Poisson Routines
Generates a single random deviate from a Poisson distribution with mean mu.

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 Floats

Implementation

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;


HTML generated by Time2HELP
http://www.time2help.com