← Back to branch summary

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

« back to all changes in this revision

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