← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to Lib/special/cdflib/cdfnbn.f

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
Import upstream version 0.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/special/cdflib/cdfnbn.f
 
1
      SUBROUTINE cdfnbn(which,p,q,s,xn,pr,ompr,status,bound)
 
2
C**********************************************************************
 
3
C
 
4
C      SUBROUTINE CDFNBN ( WHICH, P, S, XN, PR, STATUS, BOUND )
 
5
C               Cumulative Distribution Function
 
6
C               Negative BiNomial distribution
 
7
C
 
8
C
 
9
C                              Function
 
10
C
 
11
C
 
12
C     Calculates any one parameter of the negative binomial
 
13
C     distribution given values for the others.
 
14
C
 
15
C     The  cumulative  negative   binomial  distribution  returns  the
 
16
C     probability that there  will be  F or fewer failures before  the
 
17
C     XNth success in binomial trials each of which has probability of
 
18
C     success PR.
 
19
C
 
20
C     The individual term of the negative binomial is the probability of
 
21
C     S failures before XN successes and is
 
22
C          Choose( S, XN+S-1 ) * PR^(XN) * (1-PR)^S
 
23
C
 
24
C
 
25
C                              Arguments
 
26
C
 
27
C
 
28
C     WHICH --> Integer indicating which of the next four argument
 
29
C               values is to be calculated from the others.
 
30
C               Legal range: 1..4
 
31
C               iwhich = 1 : Calculate P and Q from S,XN,PR and OMPR
 
32
C               iwhich = 2 : Calculate S from P,Q,XN,PR and OMPR
 
33
C               iwhich = 3 : Calculate XN from P,Q,S,PR and OMPR
 
34
C               iwhich = 4 : Calculate PR and OMPR from P,Q,S and XN
 
35
C                    INTEGER WHICH
 
36
C
 
37
C     P <--> The cumulation from 0 to S of the  negative
 
38
C            binomial distribution.
 
39
C            Input range: [0,1].
 
40
C                    DOUBLE PRECISION P
 
41
C
 
42
C     Q <--> 1-P.
 
43
C            Input range: (0, 1].
 
44
C            P + Q = 1.0.
 
45
C                    DOUBLE PRECISION Q
 
46
C
 
47
C     S <--> The upper limit of cumulation of the binomial distribution.
 
48
C            There are F or fewer failures before the XNth success.
 
49
C            Input range: [0, +infinity).
 
50
C            Search range: [0, 1E100]
 
51
C                    DOUBLE PRECISION S
 
52
C
 
53
C     XN  <--> The number of successes.
 
54
C              Input range: [0, +infinity).
 
55
C              Search range: [0, 1E100]
 
56
C                    DOUBLE PRECISION XN
 
57
C
 
58
C     PR  <--> The probability of success in each binomial trial.
 
59
C              Input range: [0,1].
 
60
C              Search range: [0,1].
 
61
C                    DOUBLE PRECISION PR
 
62
C
 
63
C     OMPR  <--> 1-PR
 
64
C              Input range: [0,1].
 
65
C              Search range: [0,1]
 
66
C              PR + OMPR = 1.0
 
67
C                    DOUBLE PRECISION OMPR
 
68
C
 
69
C     STATUS <-- 0 if calculation completed correctly
 
70
C               -I if input parameter number I is out of range
 
71
C                1 if answer appears to be lower than lowest
 
72
C                  search bound
 
73
C                2 if answer appears to be higher than greatest
 
74
C                  search bound
 
75
C                3 if P + Q .ne. 1
 
76
C                4 if PR + OMPR .ne. 1
 
77
C                    INTEGER STATUS
 
78
C
 
79
C     BOUND <-- Undefined if STATUS is 0
 
80
C
 
81
C               Bound exceeded by parameter number I if STATUS
 
82
C               is negative.
 
83
C
 
84
C               Lower search bound if STATUS is 1.
 
85
C
 
86
C               Upper search bound if STATUS is 2.
 
87
C
 
88
C
 
89
C                              Method
 
90
C
 
91
C
 
92
C     Formula   26.5.26   of   Abramowitz  and  Stegun,  Handbook   of
 
93
C     Mathematical Functions (1966) is used  to  reduce calculation of
 
94
C     the cumulative distribution  function to that of  an  incomplete
 
95
C     beta.
 
96
C
 
97
C     Computation of other parameters involve a seach for a value that
 
98
C     produces  the desired  value  of P.   The search relies  on  the
 
99
C     monotinicity of P with the other parameter.
 
100
C
 
101
C
 
102
C**********************************************************************
 
103
C     .. Parameters ..
 
104
      DOUBLE PRECISION tol
 
105
      PARAMETER (tol=1.0D-8)
 
106
      DOUBLE PRECISION atol
 
107
      PARAMETER (atol=1.0D-50)
 
108
      DOUBLE PRECISION inf
 
109
      PARAMETER (inf=1.0D100)
 
110
      DOUBLE PRECISION one
 
111
      PARAMETER (one=1.0D0)
 
112
C     ..
 
113
C     .. Scalar Arguments ..
 
114
      DOUBLE PRECISION bound,ompr,p,pr,q,s,xn
 
115
      INTEGER status,which
 
116
C     ..
 
117
C     .. Local Scalars ..
 
118
      DOUBLE PRECISION ccum,cum,fx,pq,prompr,xhi,xlo
 
119
      LOGICAL qhi,qleft,qporq
 
120
C     ..
 
121
C     .. External Functions ..
 
122
      DOUBLE PRECISION spmpar
 
123
      EXTERNAL spmpar
 
124
C     ..
 
125
C     .. External Subroutines ..
 
126
      EXTERNAL cumnbn,dinvr,dstinv,dstzr,dzror
 
127
C     ..
 
128
C     .. Intrinsic Functions ..
 
129
      INTRINSIC abs
 
130
C     ..
 
131
      IF (.NOT. ((which.LT.1).OR. (which.GT.4))) GO TO 30
 
132
      IF (.NOT. (which.LT.1)) GO TO 10
 
133
      bound = 1.0D0
 
134
      GO TO 20
 
135
 
 
136
   10 bound = 4.0D0
 
137
   20 status = -1
 
138
      RETURN
 
139
 
 
140
   30 IF (which.EQ.1) GO TO 70
 
141
      IF (.NOT. ((p.LT.0.0D0).OR. (p.GT.1.0D0))) GO TO 60
 
142
      IF (.NOT. (p.LT.0.0D0)) GO TO 40
 
143
      bound = 0.0D0
 
144
      GO TO 50
 
145
 
 
146
   40 bound = 1.0D0
 
147
   50 status = -2
 
148
      RETURN
 
149
 
 
150
   60 CONTINUE
 
151
   70 IF (which.EQ.1) GO TO 110
 
152
      IF (.NOT. ((q.LE.0.0D0).OR. (q.GT.1.0D0))) GO TO 100
 
153
      IF (.NOT. (q.LE.0.0D0)) GO TO 80
 
154
      bound = 0.0D0
 
155
      GO TO 90
 
156
 
 
157
   80 bound = 1.0D0
 
158
   90 status = -3
 
159
      RETURN
 
160
 
 
161
  100 CONTINUE
 
162
  110 IF (which.EQ.2) GO TO 130
 
163
      IF (.NOT. (s.LT.0.0D0)) GO TO 120
 
164
      bound = 0.0D0
 
165
      status = -4
 
166
      RETURN
 
167
 
 
168
  120 CONTINUE
 
169
  130 IF (which.EQ.3) GO TO 150
 
170
      IF (.NOT. (xn.LT.0.0D0)) GO TO 140
 
171
      bound = 0.0D0
 
172
      status = -5
 
173
      RETURN
 
174
 
 
175
  140 CONTINUE
 
176
  150 IF (which.EQ.4) GO TO 190
 
177
      IF (.NOT. ((pr.LT.0.0D0).OR. (pr.GT.1.0D0))) GO TO 180
 
178
      IF (.NOT. (pr.LT.0.0D0)) GO TO 160
 
179
      bound = 0.0D0
 
180
      GO TO 170
 
181
 
 
182
  160 bound = 1.0D0
 
183
  170 status = -6
 
184
      RETURN
 
185
 
 
186
  180 CONTINUE
 
187
  190 IF (which.EQ.4) GO TO 230
 
188
      IF (.NOT. ((ompr.LT.0.0D0).OR. (ompr.GT.1.0D0))) GO TO 220
 
189
      IF (.NOT. (ompr.LT.0.0D0)) GO TO 200
 
190
      bound = 0.0D0
 
191
      GO TO 210
 
192
 
 
193
  200 bound = 1.0D0
 
194
  210 status = -7
 
195
      RETURN
 
196
 
 
197
  220 CONTINUE
 
198
  230 IF (which.EQ.1) GO TO 270
 
199
      pq = p + q
 
200
      IF (.NOT. (abs(((pq)-0.5D0)-0.5D0).GT.
 
201
     +    (3.0D0*spmpar(1)))) GO TO 260
 
202
      IF (.NOT. (pq.LT.0.0D0)) GO TO 240
 
203
      bound = 0.0D0
 
204
      GO TO 250
 
205
 
 
206
  240 bound = 1.0D0
 
207
  250 status = 3
 
208
      RETURN
 
209
 
 
210
  260 CONTINUE
 
211
  270 IF (which.EQ.4) GO TO 310
 
212
      prompr = pr + ompr
 
213
      IF (.NOT. (abs(((prompr)-0.5D0)-0.5D0).GT.
 
214
     +    (3.0D0*spmpar(1)))) GO TO 300
 
215
      IF (.NOT. (prompr.LT.0.0D0)) GO TO 280
 
216
      bound = 0.0D0
 
217
      GO TO 290
 
218
 
 
219
  280 bound = 1.0D0
 
220
  290 status = 4
 
221
      RETURN
 
222
 
 
223
  300 CONTINUE
 
224
  310 IF (.NOT. (which.EQ.1)) qporq = p .LE. q
 
225
      IF ((1).EQ. (which)) THEN
 
226
          CALL cumnbn(s,xn,pr,ompr,p,q)
 
227
          status = 0
 
228
 
 
229
      ELSE IF ((2).EQ. (which)) THEN
 
230
          s = 5.0D0
 
231
          CALL dstinv(0.0D0,inf,0.5D0,0.5D0,5.0D0,atol,tol)
 
232
          status = 0
 
233
          CALL dinvr(status,s,fx,qleft,qhi)
 
234
  320     IF (.NOT. (status.EQ.1)) GO TO 350
 
235
          CALL cumnbn(s,xn,pr,ompr,cum,ccum)
 
236
          IF (.NOT. (qporq)) GO TO 330
 
237
          fx = cum - p
 
238
          GO TO 340
 
239
 
 
240
  330     fx = ccum - q
 
241
  340     CALL dinvr(status,s,fx,qleft,qhi)
 
242
          GO TO 320
 
243
 
 
244
  350     IF (.NOT. (status.EQ.-1)) GO TO 380
 
245
          IF (.NOT. (qleft)) GO TO 360
 
246
          status = 1
 
247
          bound = 0.0D0
 
248
          GO TO 370
 
249
 
 
250
  360     status = 2
 
251
          bound = inf
 
252
  370     CONTINUE
 
253
  380     CONTINUE
 
254
 
 
255
      ELSE IF ((3).EQ. (which)) THEN
 
256
          xn = 5.0D0
 
257
          CALL dstinv(0.0D0,inf,0.5D0,0.5D0,5.0D0,atol,tol)
 
258
          status = 0
 
259
          CALL dinvr(status,xn,fx,qleft,qhi)
 
260
  390     IF (.NOT. (status.EQ.1)) GO TO 420
 
261
          CALL cumnbn(s,xn,pr,ompr,cum,ccum)
 
262
          IF (.NOT. (qporq)) GO TO 400
 
263
          fx = cum - p
 
264
          GO TO 410
 
265
 
 
266
  400     fx = ccum - q
 
267
  410     CALL dinvr(status,xn,fx,qleft,qhi)
 
268
          GO TO 390
 
269
 
 
270
  420     IF (.NOT. (status.EQ.-1)) GO TO 450
 
271
          IF (.NOT. (qleft)) GO TO 430
 
272
          status = 1
 
273
          bound = 0.0D0
 
274
          GO TO 440
 
275
 
 
276
  430     status = 2
 
277
          bound = inf
 
278
  440     CONTINUE
 
279
  450     CONTINUE
 
280
 
 
281
      ELSE IF ((4).EQ. (which)) THEN
 
282
          CALL dstzr(0.0D0,1.0D0,atol,tol)
 
283
          IF (.NOT. (qporq)) GO TO 480
 
284
          status = 0
 
285
          CALL dzror(status,pr,fx,xlo,xhi,qleft,qhi)
 
286
          ompr = one - pr
 
287
  460     IF (.NOT. (status.EQ.1)) GO TO 470
 
288
          CALL cumnbn(s,xn,pr,ompr,cum,ccum)
 
289
          fx = cum - p
 
290
          CALL dzror(status,pr,fx,xlo,xhi,qleft,qhi)
 
291
          ompr = one - pr
 
292
          GO TO 460
 
293
 
 
294
  470     GO TO 510
 
295
 
 
296
  480     status = 0
 
297
          CALL dzror(status,ompr,fx,xlo,xhi,qleft,qhi)
 
298
          pr = one - ompr
 
299
  490     IF (.NOT. (status.EQ.1)) GO TO 500
 
300
          CALL cumnbn(s,xn,pr,ompr,cum,ccum)
 
301
          fx = ccum - q
 
302
          CALL dzror(status,ompr,fx,xlo,xhi,qleft,qhi)
 
303
          pr = one - ompr
 
304
          GO TO 490
 
305
 
 
306
  500     CONTINUE
 
307
  510     IF (.NOT. (status.EQ.-1)) GO TO 540
 
308
          IF (.NOT. (qleft)) GO TO 520
 
309
          status = 1
 
310
          bound = 0.0D0
 
311
          GO TO 530
 
312
 
 
313
  520     status = 2
 
314
          bound = 1.0D0
 
315
  530     CONTINUE
 
316
  540 END IF
 
317
 
 
318
      RETURN
 
319
 
 
320
      END