| Random_Neg_Binomial Routines | 
Unit
QESBPCSRandom
| Overloaded Variants | 
| Function Random_Neg_Binomial(const sk, p: Extended): Int64; | 
| Function Random_Neg_Binomial(const sk, p: Extended; RandomGenerator: TRandomGenFunction): Int64; | 
Declaration
Function Random_Neg_Binomial(const sk, p: Extended): Int64;
Description
Adapted from Fortran 77 code from the book: Dagpunar, J. 'Principles of random variate generation' Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9
| Parameters | 
| sk | Number of Failures required (Dagpunar's words!), the "XtoY" parameter of the negative binomial. Must be positive. | 
| p | 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_Neg_Binomial (const sk, p: Extended): Int64;
begin
     Result := Random_Neg_Binomial (sk, p, DelphiRandom);
End; | 
Declaration
Function Random_Neg_Binomial(const sk, p: Extended; RandomGenerator: TRandomGenFunction): Int64;Implementation
 
| function Random_Neg_Binomial (const sk, p: Extended;
     RandomGenerator: TRandomGenFunction): Int64;
const
     h = 0.7;
var
     q, x, st, uln, v, r, s, y, g: Extended;
     k, n: Int64;
     i: Integer;
begin
     if (sk <= 0.0) or (p <= 0.0) or (p >= 1.0) then
          raise EMathError.Create (rsInvalidValue);
     q := 1.0 - p;
     x := 0.0;
     st := sk;
     if (p > h) then
     begin
          v := 1.0 / Ln (p);
          k := Trunc (st);
          for i := 1 to k do
          begin
               repeat
                    r := RandomGenerator;
               until r > 0.0;
               n := Trunc (v * Ln (r));
               x := x + n;
          end;
          st := st - k;
     end;
     s := 0.0;
     uln := -Ln (VSmall);
     if (st > -uln / Ln (q)) then
          raise EMathError.Create (rsInvalidValue);
     y := XtoY (q, st);
     g := st;
     r := RandomGenerator;
     while (y <= r) do
     begin
          r := r - y;
          s := s + 1.0;
          y := y * p * g / s;
          g := g + 1.0;
     end;
     Result := Trunc (x + s + 0.5)
End; | 
|  |