← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to routines/dcd/cdff.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/dcd/cdff.f
 
1
      SUBROUTINE cdff(which,p,q,f,dfn,dfd,status,bound)
 
2
C**********************************************************************
 
3
C
 
4
C      SUBROUTINE CDFF( WHICH, P, Q, F, DFN, DFD, STATUS, BOUND )
 
5
C               Cumulative Distribution Function
 
6
C               F distribution
 
7
C
 
8
C
 
9
C                              Function
 
10
C
 
11
C
 
12
C     Calculates any one parameter of the F distribution
 
13
C     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 F,DFN and DFD
 
23
C               iwhich = 2 : Calculate F from P,Q,DFN and DFD
 
24
C               iwhich = 3 : Calculate DFN from P,Q,F and DFD
 
25
C               iwhich = 4 : Calculate DFD from P,Q,F and DFN
 
26
C                    INTEGER WHICH
 
27
C
 
28
C       P <--> The integral from 0 to F of the f-density.
 
29
C              Input range: [0,1].
 
30
C                    DOUBLE PRECISION P
 
31
C
 
32
C       Q <--> 1-P.
 
33
C              Input range: (0, 1].
 
34
C              P + Q = 1.0.
 
35
C                    DOUBLE PRECISION Q
 
36
C
 
37
C       F <--> Upper limit of integration of the f-density.
 
38
C              Input range: [0, +infinity).
 
39
C              Search range: [0,1E300]
 
40
C                    DOUBLE PRECISION F
 
41
C
 
42
C     DFN < --> Degrees of freedom of the numerator sum of squares.
 
43
C               Input range: (0, +infinity).
 
44
C               Search range: [ 1E-300, 1E300]
 
45
C                    DOUBLE PRECISION DFN
 
46
C
 
47
C     DFD < --> Degrees of freedom of the denominator sum of squares.
 
48
C               Input range: (0, +infinity).
 
49
C               Search range: [ 1E-300, 1E300]
 
50
C                    DOUBLE PRECISION DFD
 
51
C
 
52
C     STATUS <-- 0 if calculation completed correctly
 
53
C               -I if input parameter number I is out of range
 
54
C                1 if answer appears to be lower than lowest
 
55
C                  search bound
 
56
C                2 if answer appears to be higher than greatest
 
57
C                  search bound
 
58
C                3 if P + Q .ne. 1
 
59
C                    INTEGER STATUS
 
60
C
 
61
C     BOUND <-- Undefined if STATUS is 0
 
62
C
 
63
C               Bound exceeded by parameter number I if STATUS
 
64
C               is negative.
 
65
C
 
66
C               Lower search bound if STATUS is 1.
 
67
C
 
68
C               Upper search bound if STATUS is 2.
 
69
C
 
70
C
 
71
C                              Method
 
72
C
 
73
C
 
74
C     Formula   26.6.2   of   Abramowitz   and   Stegun,  Handbook  of
 
75
C     Mathematical  Functions (1966) is used to reduce the computation
 
76
C     of the  cumulative  distribution function for the  F  variate to
 
77
C     that of an incomplete beta.
 
78
C
 
79
C     Computation of other parameters involve a seach for a value that
 
80
C     produces  the desired  value  of P.   The search relies  on  the
 
81
C     monotinicity of P with the other parameter.
 
82
C
 
83
C                              WARNING
 
84
C
 
85
C     The value of the  cumulative  F distribution is  not necessarily
 
86
C     monotone in  either degrees of freedom.  There  thus may  be two
 
87
C     values  that  provide a given CDF  value.   This routine assumes
 
88
C     monotonicity and will find an arbitrary one of the two values.
 
89
C
 
90
C**********************************************************************
 
91
C     .. Parameters ..
 
92
      DOUBLE PRECISION tol
 
93
      PARAMETER (tol=1.0D-8)
 
94
      DOUBLE PRECISION atol
 
95
      PARAMETER (atol=1.0D-50)
 
96
      DOUBLE PRECISION zero,inf
 
97
      PARAMETER (zero=1.0D-300,inf=1.0D300)
 
98
C     ..
 
99
C     .. Scalar Arguments ..
 
100
      DOUBLE PRECISION bound,dfd,dfn,f,p,q
 
101
      INTEGER status,which
 
102
C     ..
 
103
C     .. Local Scalars ..
 
104
      DOUBLE PRECISION pq,fx,cum,ccum
 
105
      LOGICAL qhi,qleft,qporq
 
106
C     ..
 
107
C     .. External Functions ..
 
108
      DOUBLE PRECISION spmpar
 
109
      EXTERNAL spmpar
 
110
C     ..
 
111
C     .. External Subroutines ..
 
112
      EXTERNAL dinvr,dstinv,cumf
 
113
C     ..
 
114
C     .. Executable Statements ..
 
115
C
 
116
C     Check arguments
 
117
C
 
118
      IF (.NOT. ((which.LT.1).OR. (which.GT.4))) GO TO 30
 
119
      IF (.NOT. (which.LT.1)) GO TO 10
 
120
      bound = 1.0D0
 
121
      GO TO 20
 
122
 
 
123
   10 bound = 4.0D0
 
124
   20 status = -1
 
125
      RETURN
 
126
 
 
127
   30 IF (which.EQ.1) GO TO 70
 
128
C
 
129
C     P
 
130
C
 
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
C
 
143
C     Q
 
144
C
 
145
      IF (.NOT. ((q.LE.0.0D0).OR. (q.GT.1.0D0))) GO TO 100
 
146
      IF (.NOT. (q.LE.0.0D0)) GO TO 80
 
147
      bound = 0.0D0
 
148
      GO TO 90
 
149
 
 
150
   80 bound = 1.0D0
 
151
   90 status = -3
 
152
      RETURN
 
153
 
 
154
  100 CONTINUE
 
155
  110 IF (which.EQ.2) GO TO 130
 
156
C
 
157
C     F
 
158
C
 
159
      IF (.NOT. (f.LT.0.0D0)) GO TO 120
 
160
      bound = 0.0D0
 
161
      status = -4
 
162
      RETURN
 
163
 
 
164
  120 CONTINUE
 
165
  130 IF (which.EQ.3) GO TO 150
 
166
C
 
167
C     DFN
 
168
C
 
169
      IF (.NOT. (dfn.LE.0.0D0)) GO TO 140
 
170
      bound = 0.0D0
 
171
      status = -5
 
172
      RETURN
 
173
 
 
174
  140 CONTINUE
 
175
  150 IF (which.EQ.4) GO TO 170
 
176
C
 
177
C     DFD
 
178
C
 
179
      IF (.NOT. (dfd.LE.0.0D0)) GO TO 160
 
180
      bound = 0.0D0
 
181
      status = -6
 
182
      RETURN
 
183
 
 
184
  160 CONTINUE
 
185
  170 IF (which.EQ.1) GO TO 210
 
186
C
 
187
C     P + Q
 
188
C
 
189
      pq = p + q
 
190
      IF (.NOT. (abs(((pq)-0.5D0)-0.5D0).GT.
 
191
     +    (3.0D0*spmpar(1)))) GO TO 200
 
192
      IF (.NOT. (pq.LT.0.0D0)) GO TO 180
 
193
      bound = 0.0D0
 
194
      GO TO 190
 
195
 
 
196
  180 bound = 1.0D0
 
197
  190 status = 3
 
198
      RETURN
 
199
 
 
200
  200 CONTINUE
 
201
  210 IF (.NOT. (which.EQ.1)) qporq = p .LE. q
 
202
C
 
203
C     Select the minimum of P or Q
 
204
C
 
205
C
 
206
C     Calculate ANSWERS
 
207
C
 
208
      IF ((1).EQ. (which)) THEN
 
209
C
 
210
C     Calculating P
 
211
C
 
212
          CALL cumf(f,dfn,dfd,p,q)
 
213
          status = 0
 
214
 
 
215
      ELSE IF ((2).EQ. (which)) THEN
 
216
C
 
217
C     Calculating F
 
218
C
 
219
          f = 5.0D0
 
220
          CALL dstinv(0.0D0,inf,0.5D0,0.5D0,5.0D0,atol,tol)
 
221
          status = 0
 
222
          CALL dinvr(status,f,fx,qleft,qhi)
 
223
  220     IF (.NOT. (status.EQ.1)) GO TO 250
 
224
          CALL cumf(f,dfn,dfd,cum,ccum)
 
225
          IF (.NOT. (qporq)) GO TO 230
 
226
          fx = cum - p
 
227
          GO TO 240
 
228
 
 
229
  230     fx = ccum - q
 
230
  240     CALL dinvr(status,f,fx,qleft,qhi)
 
231
          GO TO 220
 
232
 
 
233
  250     IF (.NOT. (status.EQ.-1)) GO TO 280
 
234
          IF (.NOT. (qleft)) GO TO 260
 
235
          status = 1
 
236
          bound = 0.0D0
 
237
          GO TO 270
 
238
 
 
239
  260     status = 2
 
240
          bound = inf
 
241
  270     CONTINUE
 
242
  280     CONTINUE
 
243
 
 
244
      ELSE IF ((3).EQ. (which)) THEN
 
245
C
 
246
C     Calculating DFN
 
247
C
 
248
          dfn = 5.0D0
 
249
          CALL dstinv(zero,inf,0.5D0,0.5D0,5.0D0,atol,tol)
 
250
          status = 0
 
251
          CALL dinvr(status,dfn,fx,qleft,qhi)
 
252
  290     IF (.NOT. (status.EQ.1)) GO TO 320
 
253
          CALL cumf(f,dfn,dfd,cum,ccum)
 
254
          IF (.NOT. (qporq)) GO TO 300
 
255
          fx = cum - p
 
256
          GO TO 310
 
257
 
 
258
  300     fx = ccum - q
 
259
  310     CALL dinvr(status,dfn,fx,qleft,qhi)
 
260
          GO TO 290
 
261
 
 
262
  320     IF (.NOT. (status.EQ.-1)) GO TO 350
 
263
          IF (.NOT. (qleft)) GO TO 330
 
264
          status = 1
 
265
          bound = zero
 
266
          GO TO 340
 
267
 
 
268
  330     status = 2
 
269
          bound = inf
 
270
  340     CONTINUE
 
271
  350     CONTINUE
 
272
 
 
273
      ELSE IF ((4).EQ. (which)) THEN
 
274
C
 
275
C     Calculating DFD
 
276
C
 
277
          dfd = 5.0D0
 
278
          CALL dstinv(zero,inf,0.5D0,0.5D0,5.0D0,atol,tol)
 
279
          status = 0
 
280
          CALL dinvr(status,dfd,fx,qleft,qhi)
 
281
  360     IF (.NOT. (status.EQ.1)) GO TO 390
 
282
          CALL cumf(f,dfn,dfd,cum,ccum)
 
283
          IF (.NOT. (qporq)) GO TO 370
 
284
          fx = cum - p
 
285
          GO TO 380
 
286
 
 
287
  370     fx = ccum - q
 
288
  380     CALL dinvr(status,dfd,fx,qleft,qhi)
 
289
          GO TO 360
 
290
 
 
291
  390     IF (.NOT. (status.EQ.-1)) GO TO 420
 
292
          IF (.NOT. (qleft)) GO TO 400
 
293
          status = 1
 
294
          bound = zero
 
295
          GO TO 410
 
296
 
 
297
  400     status = 2
 
298
          bound = inf
 
299
  410     CONTINUE
 
300
  420 END IF
 
301
 
 
302
      RETURN
 
303
 
 
304
      END