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
|
SUBROUTINE cumchn(x,df,pnonc,cum,ccum)
C***********************************************************************
C
C SUBROUTINE CUMCHN(X,DF,PNONC,CUM,CCUM)
C CUMulative of the Non-central CHi-square distribution
C
C Function
C
C Calculates the cumulative non-central chi-square
C distribution, i.e., the probability that a random variable
C which follows the non-central chi-square distribution, with
C non-centrality parameter PNONC and continuous degrees of
C freedom DF, is less than or equal to X.
C
C Arguments
C
C X --> Upper limit of integration of the non-central
C chi-square distribution.
C X is DOUBLE PRECISION
C
C DF --> Degrees of freedom of the non-central
C chi-square distribution.
C DF is DOUBLE PRECISION
C
C PNONC --> Non-centrality parameter of the non-central
C chi-square distribution.
C PNONC is DOUBLE PRECIS
C
C CUM <-- Cumulative non-central chi-square distribution.
C CUM is DOUBLE PRECISIO
C
C CCUM <-- Compliment of Cumulative non-central chi-square distribut
C CCUM is DOUBLE PRECISI
C
C
C Method
C
C Uses formula 26.4.25 of Abramowitz and Stegun, Handbook of
C Mathematical Functions, US NBS (1966) to calculate the
C non-central chi-square.
C
C Variables
C
C EPS --- Convergence criterion. The sum stops when a
C term is less than EPS*SUM.
C EPS is DOUBLE PRECISIO
C
C CCUM <-- Compliment of Cumulative non-central
C chi-square distribution.
C CCUM is DOUBLE PRECISI
C
C***********************************************************************
C
C
C .. Scalar Arguments ..
DOUBLE PRECISION ccum,cum,df,pnonc,x
C ..
C .. Local Scalars ..
DOUBLE PRECISION adj,centaj,centwt,chid2,dfd2,eps,lcntaj,lcntwt,
+ lfact,pcent,pterm,sum,sumadj,term,wt,xnonc,xx,
+ abstol
INTEGER i,icent
C ..
C .. External Functions ..
DOUBLE PRECISION alngam
EXTERNAL alngam
C ..
C .. External Subroutines ..
EXTERNAL cumchi
C ..
C .. Intrinsic Functions ..
INTRINSIC dble,exp,int,log
C ..
C .. Statement Functions ..
DOUBLE PRECISION dg
LOGICAL qsmall
C ..
C .. Data statements ..
DATA eps/1.0D-5/
DATA abstol/1.0D-300/
C ..
C .. Statement Function definitions ..
qsmall(xx) = sum .LT. abstol .OR. xx .LT. eps*sum
dg(i) = df + 2.0D0*dble(i)
C ..
C
IF (.NOT. (x.LE.0.0D0)) GO TO 10
cum = 0.0D0
ccum = 1.0D0
RETURN
10 IF (.NOT. (pnonc.LE.1.0D-10)) GO TO 20
C
C
C When non-centrality parameter is (essentially) zero,
C use cumulative chi-square distribution
C
C
CALL cumchi(x,df,cum,ccum)
RETURN
20 xnonc = pnonc/2.0D0
C***********************************************************************
C
C The following code calcualtes the weight, chi-square, and
C adjustment term for the central term in the infinite series.
C The central term is the one in which the poisson weight is
C greatest. The adjustment term is the amount that must
C be subtracted from the chi-square to move up two degrees
C of freedom.
C
C***********************************************************************
icent = int(xnonc)
IF (icent.EQ.0) icent = 1
chid2 = x/2.0D0
C
C
C Calculate central weight term
C
C
lfact = alngam(dble(icent+1))
lcntwt = -xnonc + icent*log(xnonc) - lfact
centwt = exp(lcntwt)
C
C
C Calculate central chi-square
C
C
CALL cumchi(x,dg(icent),pcent,ccum)
C
C
C Calculate central adjustment term
C
C
dfd2 = dg(icent)/2.0D0
lfact = alngam(1.0D0+dfd2)
lcntaj = dfd2*log(chid2) - chid2 - lfact
centaj = exp(lcntaj)
sum = centwt*pcent
C***********************************************************************
C
C Sum backwards from the central term towards zero.
C Quit whenever either
C (1) the zero term is reached, or
C (2) the term gets small relative to the sum, or
C
C***********************************************************************
sumadj = 0.0D0
adj = centaj
wt = centwt
i = icent
C
GO TO 40
30 IF (qsmall(term) .OR. i.EQ.0) GO TO 50
40 dfd2 = dg(i)/2.0D0
C
C
C Adjust chi-square for two fewer degrees of freedom.
C The adjusted value ends up in PTERM.
C
C
adj = adj*dfd2/chid2
sumadj = sumadj + adj
pterm = pcent + sumadj
C
C
C Adjust poisson weight for J decreased by one
C
C
wt = wt* (i/xnonc)
term = wt*pterm
sum = sum + term
i = i - 1
GO TO 30
50 sumadj = centaj
C***********************************************************************
C
C Now sum forward from the central term towards infinity.
C Quit when either
C (1) the term gets small relative to the sum, or
C
C***********************************************************************
adj = centaj
wt = centwt
i = icent
C
GO TO 70
60 IF (qsmall(term)) GO TO 80
C
C
C Update weights for next higher J
C
C
70 wt = wt* (xnonc/ (i+1))
C
C
C Calculate PTERM and add term to sum
C
C
pterm = pcent - sumadj
term = wt*pterm
sum = sum + term
C
C
C Update adjustment term for DF for next iteration
C
C
i = i + 1
dfd2 = dg(i)/2.0D0
adj = adj*chid2/dfd2
sumadj = sumadj + adj
GO TO 60
80 cum = sum
ccum = 0.5D0 + (0.5D0-cum)
C
RETURN
END
|
