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
|
SUBROUTINE cumfnc(f,dfn,dfd,pnonc,cum,ccum)
C**********************************************************************
C
C F -NON- -C-ENTRAL F DISTRIBUTION
C
C
C
C Function
C
C
C COMPUTES NONCENTRAL F DISTRIBUTION WITH DFN AND DFD
C DEGREES OF FREEDOM AND NONCENTRALITY PARAMETER PNONC
C
C
C Arguments
C
C
C X --> UPPER LIMIT OF INTEGRATION OF NONCENTRAL F IN EQUATION
C
C DFN --> DEGREES OF FREEDOM OF NUMERATOR
C
C DFD --> DEGREES OF FREEDOM OF DENOMINATOR
C
C PNONC --> NONCENTRALITY PARAMETER.
C
C CUM <-- CUMULATIVE NONCENTRAL F DISTRIBUTION
C
C CCUM <-- COMPLIMENT OF CUMMULATIVE
C
C
C Method
C
C
C USES FORMULA 26.6.20 OF REFERENCE FOR INFINITE SERIES.
C SERIES IS CALCULATED BACKWARD AND FORWARD FROM J = LAMBDA/2
C (THIS IS THE TERM WITH THE LARGEST POISSON WEIGHT) UNTIL
C THE CONVERGENCE CRITERION IS MET.
C
C FOR SPEED, THE INCOMPLETE BETA FUNCTIONS ARE EVALUATED
C BY FORMULA 26.5.16.
C
C
C REFERENCE
C
C
C HANDBOOD OF MATHEMATICAL FUNCTIONS
C EDITED BY MILTON ABRAMOWITZ AND IRENE A. STEGUN
C NATIONAL BUREAU OF STANDARDS APPLIED MATEMATICS SERIES - 55
C MARCH 1965
C P 947, EQUATIONS 26.6.17, 26.6.18
C
C
C Note
C
C
C THE SUM CONTINUES UNTIL A SUCCEEDING TERM IS LESS THAN EPS
C TIMES THE SUM (OR THE SUM IS LESS THAN 1.0E-20). EPS IS
C SET TO 1.0E-4 IN A DATA STATEMENT WHICH CAN BE CHANGED.
C
C**********************************************************************
C .. Scalar Arguments ..
DOUBLE PRECISION dfd,dfn,pnonc,f,cum,ccum
C ..
C .. Local Scalars ..
DOUBLE PRECISION dsum,dummy,prod,xx,yy
DOUBLE PRECISION adn,aup,b,betdn,betup,centwt,dnterm,eps,sum,
+ upterm,xmult,xnonc,x,abstol
INTEGER i,icent,ierr
C ..
C .. External Functions ..
DOUBLE PRECISION alngam
EXTERNAL alngam
C ..
C .. Intrinsic Functions ..
INTRINSIC log,dble,exp
C ..
C .. Statement Functions ..
LOGICAL qsmall
C ..
C .. External Subroutines ..
EXTERNAL bratio,cumf
C ..
C .. Parameters ..
DOUBLE PRECISION half
PARAMETER (half=0.5D0)
DOUBLE PRECISION done
PARAMETER (done=1.0D0)
C ..
C .. Data statements ..
DATA eps/1.0D-4/
DATA abstol/1.0D-300/
C ..
C .. Statement Function definitions ..
qsmall(x) = sum .LT. abstol .OR. x .LT. eps*sum
C ..
C .. Executable Statements ..
C
IF (.NOT. (f.LE.0.0D0)) GO TO 10
cum = 0.0D0
ccum = 1.0D0
RETURN
10 IF (.NOT. (pnonc.LT.1.0D-10)) GO TO 20
C
C Handle case in which the non-centrality parameter is
C (essentially) zero.
CALL cumf(f,dfn,dfd,cum,ccum)
RETURN
20 xnonc = pnonc/2.0D0
C Calculate the central term of the poisson weighting factor.
icent = xnonc
IF (icent.EQ.0) icent = 1
C Compute central weight term
centwt = exp(-xnonc+icent*log(xnonc)-alngam(dble(icent+1)))
C Compute central incomplete beta term
C Assure that minimum of arg to beta and 1 - arg is computed
C accurately.
prod = dfn*f
dsum = dfd + prod
yy = dfd/dsum
IF (yy.GT.half) THEN
xx = prod/dsum
yy = done - xx
ELSE
xx = done - yy
END IF
CALL bratio(dfn*half+dble(icent),dfd*half,xx,yy,betdn,dummy,ierr)
adn = dfn/2.0D0 + dble(icent)
aup = adn
b = dfd/2.0D0
betup = betdn
sum = centwt*betdn
C Now sum terms backward from icent until convergence or all done
xmult = centwt
i = icent
dnterm = exp(alngam(adn+b)-alngam(adn+1.0D0)-alngam(b)+
+ adn*log(xx)+b*log(yy))
30 IF (qsmall(xmult*betdn) .OR. i.LE.0) GO TO 40
xmult = xmult* (i/xnonc)
i = i - 1
adn = adn - 1
dnterm = (adn+1)/ ((adn+b)*xx)*dnterm
betdn = betdn + dnterm
sum = sum + xmult*betdn
GO TO 30
40 i = icent + 1
C Now sum forwards until convergence
xmult = centwt
IF ((aup-1+b).EQ.0) THEN
upterm = exp(-alngam(aup)-alngam(b)+ (aup-1)*log(xx)+
+ b*log(yy))
ELSE
upterm = exp(alngam(aup-1+b)-alngam(aup)-alngam(b)+
+ (aup-1)*log(xx)+b*log(yy))
END IF
GO TO 60
50 IF (qsmall(xmult*betup)) GO TO 70
60 xmult = xmult* (xnonc/i)
i = i + 1
aup = aup + 1
upterm = (aup+b-2.0D0)*xx/ (aup-1)*upterm
betup = betup - upterm
sum = sum + xmult*betup
GO TO 50
70 cum = sum
ccum = 0.5D0 + (0.5D0-cum)
RETURN
END
|
