1 INTEGER FUNCTION ignnbn(n,p) 2C********************************************************************** 3C 4C INTEGER FUNCTION IGNNBN( N, P ) 5C 6C GENerate Negative BiNomial random deviate 7C 8C 9C Function 10C 11C 12C Generates a single random deviate from a negative binomial 13C distribution. 14C 15C 16C Arguments 17C 18C 19C N --> Required number of events. 20C INTEGER N 21C JJV (N > 0) 22C 23C P --> The probability of an event during a Bernoulli trial. 24C REAL P 25C JJV (0.0 < P < 1.0) 26C 27C 28C 29C Method 30C 31C 32C Algorithm from page 480 of 33C 34C Devroye, Luc 35C 36C Non-Uniform Random Variate Generation. Springer-Verlag, 37C New York, 1986. 38C 39C********************************************************************** 40C .. 41C .. Scalar Arguments .. 42 REAL p 43 INTEGER n 44C .. 45C .. Local Scalars .. 46 REAL y,a,r 47C .. 48C .. External Functions .. 49C JJV changed to call SGAMMA directly 50C REAL gengam 51 REAL sgamma 52 INTEGER ignpoi 53C EXTERNAL gengam,ignpoi 54 EXTERNAL sgamma,ignpoi 55C .. 56C .. Intrinsic Functions .. 57 INTRINSIC real 58C .. 59C .. Executable Statements .. 60C Check Arguments 61C JJV changed argumnet checker to abort if N <= 0 62 IF (n.LE.0) STOP 'N <= 0 in IGNNBN' 63 IF (p.LE.0.0) STOP 'P <= 0.0 in IGNNBN' 64 IF (p.GE.1.0) STOP 'P >= 1.0 in IGNNBN' 65 66C Generate Y, a random gamma (n,(1-p)/p) variable 67C JJV Note: the above parametrization is consistent with Devroye, 68C JJV but gamma (p/(1-p),n) is the equivalent in our code 69 10 r = real(n) 70 a = p/ (1.0-p) 71C y = gengam(a,r) 72 y = sgamma(r)/a 73 74C Generate a random Poisson(y) variable 75 ignnbn = ignpoi(y) 76 RETURN 77 78 END 79