← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

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

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

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