GNU Octave  3.8.0 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
ignbin.f
Go to the documentation of this file.
1  INTEGER FUNCTION ignbin(n,pp)
2 C**********************************************************************
3 C
4 C INTEGER FUNCTION IGNBIN( N, PP )
5 C
6 C GENerate BINomial random deviate
7 C
8 C
9 C Function
10 C
11 C
12 C Generates a single random deviate from a binomial
13 C distribution whose number of trials is N and whose
14 C probability of an event in each trial is P.
15 C
16 C
17 C Arguments
18 C
19 C
20 C N --> The number of trials in the binomial distribution
21 C from which a random deviate is to be generated.
22 C INTEGER N
23 C JJV (N >= 0)
24 C
25 C PP --> The probability of an event in each trial of the
26 C binomial distribution from which a random deviate
27 C is to be generated.
28 C REAL PP
29 C JJV (0.0 <= pp <= 1.0)
30 C
31 C IGNBIN <-- A random deviate yielding the number of events
32 C from N independent trials, each of which has
33 C a probability of event P.
34 C INTEGER IGNBIN
35 C
36 C
37 C Note
38 C
39 C
40 C Uses RANF so the value of the seeds, ISEED1 and ISEED2 must be set
41 C by a call similar to the following
42 C DUM = RANSET( ISEED1, ISEED2 )
43 C
44 C
45 C Method
46 C
47 C
48 C This is algorithm BTPE from:
49 C
50 C Kachitvichyanukul, V. and Schmeiser, B. W.
51 C
52 C Binomial Random Variate Generation.
53 C Communications of the ACM, 31, 2
54 C (February, 1988) 216.
55 C
56 C**********************************************************************
57 C SUBROUTINE BTPEC(N,PP,ISEED,JX)
58 C
59 C BINOMIAL RANDOM VARIATE GENERATOR
60 C MEAN .LT. 30 -- INVERSE CDF
61 C MEAN .GE. 30 -- ALGORITHM BTPE: ACCEPTANCE-REJECTION VIA
62 C FOUR REGION COMPOSITION. THE FOUR REGIONS ARE A TRIANGLE
63 C (SYMMETRIC IN THE CENTER), A PAIR OF PARALLELOGRAMS (ABOVE
64 C THE TRIANGLE), AND EXPONENTIAL LEFT AND RIGHT TAILS.
65 C
66 C BTPE REFERS TO BINOMIAL-TRIANGLE-PARALLELOGRAM-EXPONENTIAL.
67 C BTPEC REFERS TO BTPE AND "COMBINED." THUS BTPE IS THE
68 C RESEARCH AND BTPEC IS THE IMPLEMENTATION OF A COMPLETE
69 C USABLE ALGORITHM.
70 C REFERENCE: VORATAS KACHITVICHYANUKUL AND BRUCE SCHMEISER,
71 C "BINOMIAL RANDOM VARIATE GENERATION,"
72 C COMMUNICATIONS OF THE ACM, FORTHCOMING
73 C WRITTEN: SEPTEMBER 1980.
74 C LAST REVISED: MAY 1985, JULY 1987
75 C REQUIRED SUBPROGRAM: RAND() -- A UNIFORM (0,1) RANDOM NUMBER
76 C GENERATOR
77 C ARGUMENTS
78 C
79 C N : NUMBER OF BERNOULLI TRIALS (INPUT)
80 C PP : PROBABILITY OF SUCCESS IN EACH TRIAL (INPUT)
81 C ISEED: RANDOM NUMBER SEED (INPUT AND OUTPUT)
82 C JX: RANDOMLY GENERATED OBSERVATION (OUTPUT)
83 C
84 C VARIABLES
85 C PSAVE: VALUE OF PP FROM THE LAST CALL TO BTPEC
86 C NSAVE: VALUE OF N FROM THE LAST CALL TO BTPEC
87 C XNP: VALUE OF THE MEAN FROM THE LAST CALL TO BTPEC
88 C
89 C P: PROBABILITY USED IN THE GENERATION PHASE OF BTPEC
90 C FFM: TEMPORARY VARIABLE EQUAL TO XNP + P
91 C M: INTEGER VALUE OF THE CURRENT MODE
92 C FM: FLOATING POINT VALUE OF THE CURRENT MODE
93 C XNPQ: TEMPORARY VARIABLE USED IN SETUP AND SQUEEZING STEPS
94 C P1: AREA OF THE TRIANGLE
95 C C: HEIGHT OF THE PARALLELOGRAMS
96 C XM: CENTER OF THE TRIANGLE
97 C XL: LEFT END OF THE TRIANGLE
98 C XR: RIGHT END OF THE TRIANGLE
99 C AL: TEMPORARY VARIABLE
100 C XLL: RATE FOR THE LEFT EXPONENTIAL TAIL
101 C XLR: RATE FOR THE RIGHT EXPONENTIAL TAIL
102 C P2: AREA OF THE PARALLELOGRAMS
103 C P3: AREA OF THE LEFT EXPONENTIAL TAIL
104 C P4: AREA OF THE RIGHT EXPONENTIAL TAIL
105 C U: A U(0,P4) RANDOM VARIATE USED FIRST TO SELECT ONE OF THE
106 C FOUR REGIONS AND THEN CONDITIONALLY TO GENERATE A VALUE
107 C FROM THE REGION
108 C V: A U(0,1) RANDOM NUMBER USED TO GENERATE THE RANDOM VALUE
109 C (REGION 1) OR TRANSFORMED INTO THE VARIATE TO ACCEPT OR
110 C REJECT THE CANDIDATE VALUE
111 C IX: INTEGER CANDIDATE VALUE
112 C X: PRELIMINARY CONTINUOUS CANDIDATE VALUE IN REGION 2 LOGIC
113 C AND A FLOATING POINT IX IN THE ACCEPT/REJECT LOGIC
114 C K: ABSOLUTE VALUE OF (IX-M)
115 C F: THE HEIGHT OF THE SCALED DENSITY FUNCTION USED IN THE
116 C ACCEPT/REJECT DECISION WHEN BOTH M AND IX ARE SMALL
117 C ALSO USED IN THE INVERSE TRANSFORMATION
118 C R: THE RATIO P/Q
119 C G: CONSTANT USED IN CALCULATION OF PROBABILITY
120 C MP: MODE PLUS ONE, THE LOWER INDEX FOR EXPLICIT CALCULATION
121 C OF F WHEN IX IS GREATER THAN M
122 C IX1: CANDIDATE VALUE PLUS ONE, THE LOWER INDEX FOR EXPLICIT
123 C CALCULATION OF F WHEN IX IS LESS THAN M
124 C I: INDEX FOR EXPLICIT CALCULATION OF F FOR BTPE
125 C AMAXP: MAXIMUM ERROR OF THE LOGARITHM OF NORMAL BOUND
126 C YNORM: LOGARITHM OF NORMAL BOUND
127 C ALV: NATURAL LOGARITHM OF THE ACCEPT/REJECT VARIATE V
128 C
129 C X1,F1,Z,W,Z2,X2,F2, AND W2 ARE TEMPORARY VARIABLES TO BE
130 C USED IN THE FINAL ACCEPT/REJECT TEST
131 C
132 C QN: PROBABILITY OF NO SUCCESS IN N TRIALS
133 C
134 C REMARK
135 C IX AND JX COULD LOGICALLY BE THE SAME VARIABLE, WHICH WOULD
136 C SAVE A MEMORY POSITION AND A LINE OF CODE. HOWEVER, SOME
137 C COMPILERS (E.G.,CDC MNF) OPTIMIZE BETTER WHEN THE ARGUMENTS
138 C ARE NOT INVOLVED.
139 C
140 C ISEED NEEDS TO BE DOUBLE PRECISION IF THE IMSL ROUTINE
141 C GGUBFS IS USED TO GENERATE UNIFORM RANDOM NUMBER, OTHERWISE
142 C TYPE OF ISEED SHOULD BE DICTATED BY THE UNIFORM GENERATOR
143 C
144 C**********************************************************************
145
146 C
147 C
148 C
149 C*****DETERMINE APPROPRIATE ALGORITHM AND WHETHER SETUP IS NECESSARY
150 C
151 C ..
152 C .. Scalar Arguments ..
153  REAL pp
154  INTEGER n
155 C ..
156 C .. Local Scalars ..
157  REAL al,alv,amaxp,c,f,f1,f2,ffm,fm,g,p,p1,p2,p3,p4,psave,q,qn,r,u,
158  + v,w,w2,x,x1,x2,xl,xll,xlr,xm,xnp,xnpq,xr,ynorm,z,z2
159  INTEGER i,ix,ix1,k,m,mp,nsave
160 C ..
161 C .. External Functions ..
162  REAL ranf
163  EXTERNAL ranf
164 C ..
165 C .. Intrinsic Functions ..
166  INTRINSIC abs,alog,amin1,iabs,int,sqrt
167 C JJV ..
168 C JJV .. Save statement ..
169  SAVE p,q,m,fm,xnp,xnpq,p1,xm,xl,xr,c,xll,xlr,p2,p3,p4,qn,r,g,
170  + psave,nsave
171 C JJV I am including the variables in data statements
172 C ..
173 C .. Data statements ..
174 C JJV made these ridiculous starting values - the hope is that
175 C JJV no one will call this the first time with them as args
176  DATA psave,nsave/-1.0e37,-214748365/
177 C ..
178 C .. Executable Statements ..
179  IF (pp.NE.psave) go to 10
180  IF (n.NE.nsave) go to 20
181  IF (xnp-30.0.LT.0.0) go to 150
182  go to 30
183 C
184 C*****SETUP, PERFORM ONLY WHEN PARAMETERS CHANGE
185 C
186
187 C JJV added the argument checker - involved only renaming 10
188 C JJV and 20 to the checkers and adding checkers
189 C JJV Only remaining problem - if called initially with the
190 C JJV initial values of psave and nsave, it will hang
191  10 IF (pp.LT.0.0) CALL xstopx('PP < 0.0 in IGNBIN - ABORT!')
192  IF (pp.GT.1.0) CALL xstopx('PP > 1.0 in IGNBIN - ABORT!')
193  psave = pp
194  p = amin1(psave,1.-psave)
195  q = 1. - p
196  20 IF (n.LT.0) CALL xstopx('N < 0 in IGNBIN - ABORT!')
197  xnp = n*p
198  nsave = n
199  IF (xnp.LT.30.) go to 140
200  ffm = xnp + p
201  m = ffm
202  fm = m
203  xnpq = xnp*q
204  p1 = int(2.195*sqrt(xnpq)-4.6*q) + 0.5
205  xm = fm + 0.5
206  xl = xm - p1
207  xr = xm + p1
208  c = 0.134 + 20.5/ (15.3+fm)
209  al = (ffm-xl)/ (ffm-xl*p)
210  xll = al* (1.+.5*al)
211  al = (xr-ffm)/ (xr*q)
212  xlr = al* (1.+.5*al)
213  p2 = p1* (1.+c+c)
214  p3 = p2 + c/xll
215  p4 = p3 + c/xlr
216 C WRITE(6,100) N,P,P1,P2,P3,P4,XL,XR,XM,FM
217 C 100 FORMAT(I15,4F18.7/5F18.7)
218 C
219 C*****GENERATE VARIATE
220 C
221  30 u = ranf()*p4
222  v = ranf()
223 C
224 C TRIANGULAR REGION
225 C
226  IF (u.GT.p1) go to 40
227  ix = xm - p1*v + u
228  go to 170
229 C
230 C PARALLELOGRAM REGION
231 C
232  40 IF (u.GT.p2) go to 50
233  x = xl + (u-p1)/c
234  v = v*c + 1. - abs(xm-x)/p1
235  IF (v.GT.1. .OR. v.LE.0.) go to 30
236  ix = x
237  go to 70
238 C
239 C LEFT TAIL
240 C
241  50 IF (u.GT.p3) go to 60
242  ix = xl + alog(v)/xll
243  IF (ix.LT.0) go to 30
244  v = v* (u-p2)*xll
245  go to 70
246 C
247 C RIGHT TAIL
248 C
249  60 ix = xr - alog(v)/xlr
250  IF (ix.GT.n) go to 30
251  v = v* (u-p3)*xlr
252 C
253 C*****DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST
254 C
255  70 k = iabs(ix-m)
256  IF (k.GT.20 .AND. k.LT.xnpq/2-1) go to 130
257 C
258 C EXPLICIT EVALUATION
259 C
260  f = 1.0
261  r = p/q
262  g = (n+1)*r
263  IF (m-ix.LT.0) go to 80
264  IF (m-ix.EQ.0) go to 120
265  go to 100
266  80 mp = m + 1
267  DO 90 i = mp,ix
268  f = f* (g/i-r)
269  90 CONTINUE
270  go to 120
271
272  100 ix1 = ix + 1
273  DO 110 i = ix1,m
274  f = f/ (g/i-r)
275  110 CONTINUE
276  120 IF (v-f.LE.0) go to 170
277  go to 30
278 C
279 C SQUEEZING USING UPPER AND LOWER BOUNDS ON ALOG(F(X))
280 C
281  130 amaxp = (k/xnpq)* ((k* (k/3.+.625)+.1666666666666)/xnpq+.5)
282  ynorm = -k*k/ (2.*xnpq)
283  alv = alog(v)
284  IF (alv.LT.ynorm-amaxp) go to 170
285  IF (alv.GT.ynorm+amaxp) go to 30
286 C
287 C STIRLING'S FORMULA TO MACHINE ACCURACY FOR
288 C THE FINAL ACCEPTANCE/REJECTION TEST
289 C
290  x1 = ix + 1
291  f1 = fm + 1.
292  z = n + 1 - fm
293  w = n - ix + 1.
294  z2 = z*z
295  x2 = x1*x1
296  f2 = f1*f1
297  w2 = w*w
298  IF (alv- (xm*alog(f1/x1)+ (n-m+.5)*alog(z/w)+ (ix-
299  + m)*alog(w*p/ (x1*q))+ (13860.- (462.- (132.- (99.-
300  + 140./f2)/f2)/f2)/f2)/f1/166320.+ (13860.- (462.- (132.- (99.-
301  + 140./z2)/z2)/z2)/z2)/z/166320.+ (13860.- (462.- (132.- (99.-
302  + 140./x2)/x2)/x2)/x2)/x1/166320.+ (13860.- (462.- (132.- (99.-
303  + 140./w2)/w2)/w2)/w2)/w/166320.) .LE. 0.) go to 170
304  go to 30
305 C
306 C INVERSE CDF LOGIC FOR MEAN LESS THAN 30
307 C
308  140 qn = q**n
309  r = p/q
310  g = r* (n+1)
311  150 ix = 0
312  f = qn
313  u = ranf()
314  160 IF (u.LT.f) go to 170
315  IF (ix.GT.110) go to 150
316  u = u - f
317  ix = ix + 1
318  f = f* (g/ix-r)
319  go to 160
320
321  170 IF (psave.GT.0.5) ix = n - ix
322  ignbin = ix
323  RETURN
324
325  END