GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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
Generated on Mon Dec 30 2013 03:04:47 for GNU Octave by   doxygen 1.8.1.2