File: cdff.f

package info (click to toggle)
python-scipy 0.18.1-2
  • links: PTS, VCS
  • area: main
  • in suites: stretch
  • size: 75,464 kB
  • ctags: 79,406
  • sloc: python: 143,495; cpp: 89,357; fortran: 81,650; ansic: 79,778; makefile: 364; sh: 265
file content (267 lines) | stat: -rw-r--r-- 7,465 bytes parent folder | download | duplicates (8)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
      SUBROUTINE cdff(which,p,q,f,dfn,dfd,status,bound)
C**********************************************************************
C
C      SUBROUTINE CDFF( WHICH, P, Q, F, DFN, DFD, STATUS, BOUND )
C               Cumulative Distribution Function
C               F distribution
C
C
C                              Function
C
C
C     Calculates any one parameter of the F distribution
C     given values for the others.
C
C
C                              Arguments
C
C
C     WHICH --> Integer indicating which of the next four argument
C               values is to be calculated from the others.
C               Legal range: 1..4
C               iwhich = 1 : Calculate P and Q from F,DFN and DFD
C               iwhich = 2 : Calculate F from P,Q,DFN and DFD
C               iwhich = 3 : Calculate DFN from P,Q,F and DFD
C               iwhich = 4 : Calculate DFD from P,Q,F and DFN
C                    INTEGER WHICH
C
C       P <--> The integral from 0 to F of the f-density.
C              Input range: [0,1].
C                    DOUBLE PRECISION P
C
C       Q <--> 1-P.
C              Input range: (0, 1].
C              P + Q = 1.0.
C                    DOUBLE PRECISION Q
C
C       F <--> Upper limit of integration of the f-density.
C              Input range: [0, +infinity).
C              Search range: [0,1E100]
C                    DOUBLE PRECISION F
C
C     DFN < --> Degrees of freedom of the numerator sum of squares.
C               Input range: (0, +infinity).
C               Search range: [ 1E-100, 1E100]
C                    DOUBLE PRECISION DFN
C
C     DFD < --> Degrees of freedom of the denominator sum of squares.
C               Input range: (0, +infinity).
C               Search range: [ 1E-100, 1E100]
C                    DOUBLE PRECISION DFD
C
C     STATUS <-- 0 if calculation completed correctly
C               -I if input parameter number I is out of range
C                1 if answer appears to be lower than lowest
C                  search bound
C                2 if answer appears to be higher than greatest
C                  search bound
C                3 if P + Q .ne. 1
C                    INTEGER STATUS
C
C     BOUND <-- Undefined if STATUS is 0
C
C               Bound exceeded by parameter number I if STATUS
C               is negative.
C
C               Lower search bound if STATUS is 1.
C
C               Upper search bound if STATUS is 2.
C
C
C                              Method
C
C
C     Formula   26.6.2   of   Abramowitz   and   Stegun,  Handbook  of
C     Mathematical  Functions (1966) is used to reduce the computation
C     of the  cumulative  distribution function for the  F  variate to
C     that of an incomplete beta.
C
C     Computation of other parameters involve a seach for a value that
C     produces  the desired  value  of P.   The search relies  on  the
C     monotinicity of P with the other parameter.
C
C                              WARNING
C
C     The value of the  cumulative  F distribution is  not necessarily
C     monotone in  either degrees of freedom.  There  thus may  be two
C     values  that  provide a given CDF  value.   This routine assumes
C     monotonicity and will find an arbitrary one of the two values.
C
C**********************************************************************
C     .. Parameters ..
      DOUBLE PRECISION tol
      PARAMETER (tol=1.0D-8)
      DOUBLE PRECISION atol
      PARAMETER (atol=1.0D-50)
      DOUBLE PRECISION zero,inf
      PARAMETER (zero=1.0D-100,inf=1.0D100)
C     ..
C     .. Scalar Arguments ..
      DOUBLE PRECISION bound,dfd,dfn,f,p,q
      INTEGER status,which
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION ccum,cum,fx,pq
      LOGICAL qhi,qleft,qporq
C     ..
C     .. External Functions ..
      DOUBLE PRECISION spmpar
      EXTERNAL spmpar
C     ..
C     .. External Subroutines ..
      EXTERNAL cumf,dinvr,dstinv
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC abs
C     ..
      IF (.NOT. ((which.LT.1).OR. (which.GT.4))) GO TO 30
      IF (.NOT. (which.LT.1)) GO TO 10
      bound = 1.0D0
      GO TO 20

   10 bound = 4.0D0
   20 status = -1
      RETURN

   30 IF (which.EQ.1) GO TO 70
      IF (.NOT. ((p.LT.0.0D0).OR. (p.GT.1.0D0))) GO TO 60
      IF (.NOT. (p.LT.0.0D0)) GO TO 40
      bound = 0.0D0
      GO TO 50

   40 bound = 1.0D0
   50 status = -2
      RETURN

   60 CONTINUE
   70 IF (which.EQ.1) GO TO 110
      IF (.NOT. ((q.LE.0.0D0).OR. (q.GT.1.0D0))) GO TO 100
      IF (.NOT. (q.LE.0.0D0)) GO TO 80
      bound = 0.0D0
      GO TO 90

   80 bound = 1.0D0
   90 status = -3
      RETURN

  100 CONTINUE
  110 IF (which.EQ.2) GO TO 130
      IF (.NOT. (f.LT.0.0D0)) GO TO 120
      bound = 0.0D0
      status = -4
      RETURN

  120 CONTINUE
  130 IF (which.EQ.3) GO TO 150
      IF (.NOT. (dfn.LE.0.0D0)) GO TO 140
      bound = 0.0D0
      status = -5
      RETURN

  140 CONTINUE
  150 IF (which.EQ.4) GO TO 170
      IF (.NOT. (dfd.LE.0.0D0)) GO TO 160
      bound = 0.0D0
      status = -6
      RETURN

  160 CONTINUE
  170 IF (which.EQ.1) GO TO 210
      pq = p + q
      IF (.NOT. (abs(((pq)-0.5D0)-0.5D0).GT.
     +    (3.0D0*spmpar(1)))) GO TO 200
      IF (.NOT. (pq.LT.0.0D0)) GO TO 180
      bound = 0.0D0
      GO TO 190

  180 bound = 1.0D0
  190 status = 3
      RETURN

  200 CONTINUE
  210 IF (.NOT. (which.EQ.1)) qporq = p .LE. q
      IF ((1).EQ. (which)) THEN
          CALL cumf(f,dfn,dfd,p,q)
          status = 0

      ELSE IF ((2).EQ. (which)) THEN
          f = 5.0D0
          CALL dstinv(0.0D0,inf,0.5D0,0.5D0,5.0D0,atol,tol)
          status = 0
          CALL dinvr(status,f,fx,qleft,qhi)
  220     IF (.NOT. (status.EQ.1)) GO TO 250
          CALL cumf(f,dfn,dfd,cum,ccum)
          IF (.NOT. (qporq)) GO TO 230
          fx = cum - p
          GO TO 240

  230     fx = ccum - q
  240     CALL dinvr(status,f,fx,qleft,qhi)
          GO TO 220

  250     IF (.NOT. (status.EQ.-1)) GO TO 280
          IF (.NOT. (qleft)) GO TO 260
          status = 1
          bound = 0.0D0
          GO TO 270

  260     status = 2
          bound = inf
  270     CONTINUE
  280     CONTINUE

      ELSE IF ((3).EQ. (which)) THEN
          dfn = 5.0D0
          CALL dstinv(zero,inf,0.5D0,0.5D0,5.0D0,atol,tol)
          status = 0
          CALL dinvr(status,dfn,fx,qleft,qhi)
  290     IF (.NOT. (status.EQ.1)) GO TO 320
          CALL cumf(f,dfn,dfd,cum,ccum)
          IF (.NOT. (qporq)) GO TO 300
          fx = cum - p
          GO TO 310

  300     fx = ccum - q
  310     CALL dinvr(status,dfn,fx,qleft,qhi)
          GO TO 290

  320     IF (.NOT. (status.EQ.-1)) GO TO 350
          IF (.NOT. (qleft)) GO TO 330
          status = 1
          bound = zero
          GO TO 340

  330     status = 2
          bound = inf
  340     CONTINUE
  350     CONTINUE

      ELSE IF ((4).EQ. (which)) THEN
          dfd = 5.0D0
          CALL dstinv(zero,inf,0.5D0,0.5D0,5.0D0,atol,tol)
          status = 0
          CALL dinvr(status,dfd,fx,qleft,qhi)
  360     IF (.NOT. (status.EQ.1)) GO TO 390
          CALL cumf(f,dfn,dfd,cum,ccum)
          IF (.NOT. (qporq)) GO TO 370
          fx = cum - p
          GO TO 380

  370     fx = ccum - q
  380     CALL dinvr(status,dfd,fx,qleft,qhi)
          GO TO 360

  390     IF (.NOT. (status.EQ.-1)) GO TO 420
          IF (.NOT. (qleft)) GO TO 400
          status = 1
          bound = zero
          GO TO 410

  400     status = 2
          bound = inf
  410     CONTINUE
  420 END IF

      RETURN

      END