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
|