00001 INTEGER FUNCTION ignnbn(n,p) 00002 C********************************************************************** 00003 C 00004 C INTEGER FUNCTION IGNNBN( N, P ) 00005 C 00006 C GENerate Negative BiNomial random deviate 00007 C 00008 C 00009 C Function 00010 C 00011 C 00012 C Generates a single random deviate from a negative binomial 00013 C distribution. 00014 C 00015 C 00016 C Arguments 00017 C 00018 C 00019 C N --> Required number of events. 00020 C INTEGER N 00021 C JJV (N > 0) 00022 C 00023 C P --> The probability of an event during a Bernoulli trial. 00024 C REAL P 00025 C JJV (0.0 < P < 1.0) 00026 C 00027 C 00028 C 00029 C Method 00030 C 00031 C 00032 C Algorithm from page 480 of 00033 C 00034 C Devroye, Luc 00035 C 00036 C Non-Uniform Random Variate Generation. Springer-Verlag, 00037 C New York, 1986. 00038 C 00039 C********************************************************************** 00040 C .. 00041 C .. Scalar Arguments .. 00042 REAL p 00043 INTEGER n 00044 C .. 00045 C .. Local Scalars .. 00046 REAL y,a,r 00047 C .. 00048 C .. External Functions .. 00049 C JJV changed to call SGAMMA directly 00050 C REAL gengam 00051 REAL sgamma 00052 INTEGER ignpoi 00053 C EXTERNAL gengam,ignpoi 00054 EXTERNAL sgamma,ignpoi 00055 C .. 00056 C .. Intrinsic Functions .. 00057 INTRINSIC real 00058 C .. 00059 C .. Executable Statements .. 00060 C Check Arguments 00061 C JJV changed argumnet checker to abort if N <= 0 00062 IF (n.LE.0) CALL XSTOPX ('N <= 0 in IGNNBN') 00063 IF (p.LE.0.0) CALL XSTOPX ('P <= 0.0 in IGNNBN') 00064 IF (p.GE.1.0) CALL XSTOPX ('P >= 1.0 in IGNNBN') 00065 00066 C Generate Y, a random gamma (n,(1-p)/p) variable 00067 C JJV Note: the above parametrization is consistent with Devroye, 00068 C JJV but gamma (p/(1-p),n) is the equivalent in our code 00069 10 r = real(n) 00070 a = p/ (1.0-p) 00071 C y = gengam(a,r) 00072 y = sgamma(r)/a 00073 00074 C Generate a random Poisson(y) variable 00075 ignnbn = ignpoi(y) 00076 RETURN 00077 00078 END
1.7.1
