← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/special/cdflib/cdffnc.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/cdffnc.f
 
1
      SUBROUTINE cdffnc(which,p,q,f,dfn,dfd,phonc,status,bound)
 
2
C**********************************************************************
 
3
C
 
4
C      SUBROUTINE CDFFNC( WHICH, P, Q, F, DFN, DFD, PNONC, STATUS, BOUND )
 
5
C               Cumulative Distribution Function
 
6
C               Non-central F distribution
 
7
C
 
8
C
 
9
C                              Function
 
10
C
 
11
C
 
12
C     Calculates any one parameter of the Non-central F
 
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 five argument
 
20
C               values is to be calculated from the others.
 
21
C               Legal range: 1..5
 
22
C               iwhich = 1 : Calculate P and Q from F,DFN,DFD and PNONC
 
23
C               iwhich = 2 : Calculate F from P,Q,DFN,DFD and PNONC
 
24
C               iwhich = 3 : Calculate DFN from P,Q,F,DFD and PNONC
 
25
C               iwhich = 4 : Calculate DFD from P,Q,F,DFN and PNONC
 
26
C               iwhich = 5 : Calculate PNONC from P,Q,F,DFN and DFD
 
27
C                    INTEGER WHICH
 
28
C
 
29
C       P <--> The integral from 0 to F of the non-central f-density.
 
30
C              Input range: [0,1-1E-16).
 
31
C                    DOUBLE PRECISION P
 
32
C
 
33
C       Q <--> 1-P.
 
34
C            Q is not used by this subroutine and is only included
 
35
C            for similarity with other cdf* routines.
 
36
C                    DOUBLE PRECISION Q
 
37
C
 
38
C       F <--> Upper limit of integration of the non-central f-density.
 
39
C              Input range: [0, +infinity).
 
40
C              Search range: [0,1E100]
 
41
C                    DOUBLE PRECISION F
 
42
C
 
43
C     DFN < --> Degrees of freedom of the numerator sum of squares.
 
44
C               Input range: (0, +infinity).
 
45
C               Search range: [ 1E-100, 1E100]
 
46
C                    DOUBLE PRECISION DFN
 
47
C
 
48
C     DFD < --> Degrees of freedom of the denominator sum of squares.
 
49
C               Must be in range: (0, +infinity).
 
50
C               Input range: (0, +infinity).
 
51
C               Search range: [ 1E-100, 1E100]
 
52
C                    DOUBLE PRECISION DFD
 
53
C
 
54
C     PNONC <-> The non-centrality parameter
 
55
C               Input range: [0,infinity)
 
56
C               Search range: [0,1E4]
 
57
C                    DOUBLE PRECISION PHONC
 
58
C
 
59
C     STATUS <-- 0 if calculation completed correctly
 
60
C               -I if input parameter number I is out of range
 
61
C                1 if answer appears to be lower than lowest
 
62
C                  search bound
 
63
C                2 if answer appears to be higher than greatest
 
64
C                  search bound
 
65
C                3 if P + Q .ne. 1
 
66
C                    INTEGER STATUS
 
67
C
 
68
C     BOUND <-- Undefined if STATUS is 0
 
69
C
 
70
C               Bound exceeded by parameter number I if STATUS
 
71
C               is negative.
 
72
C
 
73
C               Lower search bound if STATUS is 1.
 
74
C
 
75
C               Upper search bound if STATUS is 2.
 
76
C
 
77
C
 
78
C                              Method
 
79
C
 
80
C
 
81
C     Formula  26.6.20   of   Abramowitz   and   Stegun,  Handbook  of
 
82
C     Mathematical  Functions (1966) is used to compute the cumulative
 
83
C     distribution function.
 
84
C
 
85
C     Computation of other parameters involve a seach for a value that
 
86
C     produces  the desired  value  of P.   The search relies  on  the
 
87
C     monotinicity of P with the other parameter.
 
88
C
 
89
C                            WARNING
 
90
C
 
91
C     The computation time  required for this  routine is proportional
 
92
C     to the noncentrality  parameter  (PNONC).  Very large  values of
 
93
C     this parameter can consume immense  computer resources.  This is
 
94
C     why the search range is bounded by 10,000.
 
95
C
 
96
C                              WARNING
 
97
C
 
98
C     The  value  of the  cumulative  noncentral F distribution is not
 
99
C     necessarily monotone in either degrees  of freedom.  There  thus
 
100
C     may be two values that provide a given  CDF value.  This routine
 
101
C     assumes monotonicity  and will find  an arbitrary one of the two
 
102
C     values.
 
103
C
 
104
C**********************************************************************
 
105
C     .. Parameters ..
 
106
      DOUBLE PRECISION tent4
 
107
      PARAMETER (tent4=1.0D4)
 
108
      DOUBLE PRECISION tol
 
109
      PARAMETER (tol=1.0D-8)
 
110
      DOUBLE PRECISION atol
 
111
      PARAMETER (atol=1.0D-50)
 
112
      DOUBLE PRECISION zero,one,inf
 
113
      PARAMETER (zero=1.0D-100,one=1.0D0-1.0D-16,inf=1.0D100)
 
114
C     ..
 
115
C     .. Scalar Arguments ..
 
116
      DOUBLE PRECISION bound,dfd,dfn,f,p,phonc,q
 
117
      INTEGER status,which
 
118
C     ..
 
119
C     .. Local Scalars ..
 
120
      DOUBLE PRECISION ccum,cum,fx
 
121
      LOGICAL qhi,qleft
 
122
C     ..
 
123
C     .. External Subroutines ..
 
124
      EXTERNAL cumfnc,dinvr,dstinv
 
125
C     ..
 
126
      IF (.NOT. ((which.LT.1).OR. (which.GT.5))) GO TO 30
 
127
      IF (.NOT. (which.LT.1)) GO TO 10
 
128
      bound = 1.0D0
 
129
      GO TO 20
 
130
 
 
131
   10 bound = 5.0D0
 
132
   20 status = -1
 
133
      RETURN
 
134
 
 
135
   30 IF (which.EQ.1) GO TO 70
 
136
      IF (.NOT. ((p.LT.0.0D0).OR. (p.GT.one))) GO TO 60
 
137
      IF (.NOT. (p.LT.0.0D0)) GO TO 40
 
138
      bound = 0.0D0
 
139
      GO TO 50
 
140
 
 
141
   40 bound = one
 
142
   50 status = -2
 
143
      RETURN
 
144
 
 
145
   60 CONTINUE
 
146
   70 IF (which.EQ.2) GO TO 90
 
147
      IF (.NOT. (f.LT.0.0D0)) GO TO 80
 
148
      bound = 0.0D0
 
149
      status = -4
 
150
      RETURN
 
151
 
 
152
   80 CONTINUE
 
153
   90 IF (which.EQ.3) GO TO 110
 
154
      IF (.NOT. (dfn.LE.0.0D0)) GO TO 100
 
155
      bound = 0.0D0
 
156
      status = -5
 
157
      RETURN
 
158
 
 
159
  100 CONTINUE
 
160
  110 IF (which.EQ.4) GO TO 130
 
161
      IF (.NOT. (dfd.LE.0.0D0)) GO TO 120
 
162
      bound = 0.0D0
 
163
      status = -6
 
164
      RETURN
 
165
 
 
166
  120 CONTINUE
 
167
  130 IF (which.EQ.5) GO TO 150
 
168
      IF (.NOT. (phonc.LT.0.0D0)) GO TO 140
 
169
      bound = 0.0D0
 
170
      status = -7
 
171
      RETURN
 
172
 
 
173
  140 CONTINUE
 
174
  150 IF ((1).EQ. (which)) THEN
 
175
          CALL cumfnc(f,dfn,dfd,phonc,p,q)
 
176
          status = 0
 
177
 
 
178
      ELSE IF ((2).EQ. (which)) THEN
 
179
          f = 5.0D0
 
180
          CALL dstinv(0.0D0,inf,0.5D0,0.5D0,5.0D0,atol,tol)
 
181
          status = 0
 
182
          CALL dinvr(status,f,fx,qleft,qhi)
 
183
  160     IF (.NOT. (status.EQ.1)) GO TO 170
 
184
          CALL cumfnc(f,dfn,dfd,phonc,cum,ccum)
 
185
          fx = cum - p
 
186
          CALL dinvr(status,f,fx,qleft,qhi)
 
187
          GO TO 160
 
188
 
 
189
  170     IF (.NOT. (status.EQ.-1)) GO TO 200
 
190
          IF (.NOT. (qleft)) GO TO 180
 
191
          status = 1
 
192
          bound = 0.0D0
 
193
          GO TO 190
 
194
 
 
195
  180     status = 2
 
196
          bound = inf
 
197
  190     CONTINUE
 
198
  200     CONTINUE
 
199
 
 
200
      ELSE IF ((3).EQ. (which)) THEN
 
201
          dfn = 5.0D0
 
202
          CALL dstinv(zero,inf,0.5D0,0.5D0,5.0D0,atol,tol)
 
203
          status = 0
 
204
          CALL dinvr(status,dfn,fx,qleft,qhi)
 
205
  210     IF (.NOT. (status.EQ.1)) GO TO 220
 
206
          CALL cumfnc(f,dfn,dfd,phonc,cum,ccum)
 
207
          fx = cum - p
 
208
          CALL dinvr(status,dfn,fx,qleft,qhi)
 
209
          GO TO 210
 
210
 
 
211
  220     IF (.NOT. (status.EQ.-1)) GO TO 250
 
212
          IF (.NOT. (qleft)) GO TO 230
 
213
          status = 1
 
214
          bound = zero
 
215
          GO TO 240
 
216
 
 
217
  230     status = 2
 
218
          bound = inf
 
219
  240     CONTINUE
 
220
  250     CONTINUE
 
221
 
 
222
      ELSE IF ((4).EQ. (which)) THEN
 
223
          dfd = 5.0D0
 
224
          CALL dstinv(zero,inf,0.5D0,0.5D0,5.0D0,atol,tol)
 
225
          status = 0
 
226
          CALL dinvr(status,dfd,fx,qleft,qhi)
 
227
  260     IF (.NOT. (status.EQ.1)) GO TO 270
 
228
          CALL cumfnc(f,dfn,dfd,phonc,cum,ccum)
 
229
          fx = cum - p
 
230
          CALL dinvr(status,dfd,fx,qleft,qhi)
 
231
          GO TO 260
 
232
 
 
233
  270     IF (.NOT. (status.EQ.-1)) GO TO 300
 
234
          IF (.NOT. (qleft)) GO TO 280
 
235
          status = 1
 
236
          bound = zero
 
237
          GO TO 290
 
238
 
 
239
  280     status = 2
 
240
          bound = inf
 
241
  290     CONTINUE
 
242
  300     CONTINUE
 
243
 
 
244
      ELSE IF ((5).EQ. (which)) THEN
 
245
          phonc = 5.0D0
 
246
          CALL dstinv(0.0D0,tent4,0.5D0,0.5D0,5.0D0,atol,tol)
 
247
          status = 0
 
248
          CALL dinvr(status,phonc,fx,qleft,qhi)
 
249
  310     IF (.NOT. (status.EQ.1)) GO TO 320
 
250
          CALL cumfnc(f,dfn,dfd,phonc,cum,ccum)
 
251
          fx = cum - p
 
252
          CALL dinvr(status,phonc,fx,qleft,qhi)
 
253
          GO TO 310
 
254
 
 
255
  320     IF (.NOT. (status.EQ.-1)) GO TO 350
 
256
          IF (.NOT. (qleft)) GO TO 330
 
257
          status = 1
 
258
          bound = 0.0D0
 
259
          GO TO 340
 
260
 
 
261
  330     status = 2
 
262
          bound = tent4
 
263
  340     CONTINUE
 
264
  350 END IF
 
265
 
 
266
      RETURN
 
267
 
 
268
      END