1       SUBROUTINE cdfgam(which,p,q,x,shape,scale,status,bound)
2 C**********************************************************************
3 C
4 C      SUBROUTINE CDFGAM( WHICH, P, Q, X, SHAPE, SCALE, STATUS, BOUND )
5 C               Cumulative Distribution Function
6 C                         GAMma Distribution
7 C
8 C
9 C                              Function
10 C
11 C
12 C     Calculates any one parameter of the gamma
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 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 X,SHAPE and SCALE
23 C               iwhich = 2 : Calculate X from P,Q,SHAPE and SCALE
24 C               iwhich = 3 : Calculate SHAPE from P,Q,X and SCALE
25 C               iwhich = 4 : Calculate SCALE from P,Q,X and SHAPE
26 C                    INTEGER WHICH
27 C
28 C     P <--> The integral from 0 to X of the gamma 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
38 C     X <--> The upper limit of integration of the gamma density.
39 C            Input range: [0, +infinity).
40 C            Search range: [0,1E300]
41 C                    DOUBLE PRECISION X
42 C
43 C     SHAPE <--> The shape parameter of the gamma density.
44 C                Input range: (0, +infinity).
45 C                Search range: [1E-300,1E300]
46 C                  DOUBLE PRECISION SHAPE
47 C
48 C
49 C     SCALE <--> The scale parameter of the gamma density.
50 C                Input range: (0, +infinity).
51 C                Search range: (1E-300,1E300]
52 C                   DOUBLE PRECISION SCALE
53 C
54 C     STATUS <-- 0 if calculation completed correctly
55 C               -I if input parameter number I is out of range
56 C                1 if answer appears to be lower than lowest
57 C                  search bound
58 C                2 if answer appears to be higher than greatest
59 C                  search bound
60 C                3 if P + Q .ne. 1
61 C                10 if the gamma or inverse gamma routine cannot
62 C                   compute the answer.  Usually happens only for
63 C                   X and SHAPE very large (gt 1E10 or more)
64 C                    INTEGER STATUS
65 C
66 C     BOUND <-- Undefined if STATUS is 0
67 C
68 C               Bound exceeded by parameter number I if STATUS
69 C               is negative.
70 C
71 C               Lower search bound if STATUS is 1.
72 C
73 C               Upper search bound if STATUS is 2.
74 C
75 C
76 C                              Method
77 C
78 C
79 C     Cumulative distribution function (P) is calculated directly by
80 C     the code associated with:
81 C
82 C     DiDinato, A. R. and Morris, A. H. Computation of the  incomplete
83 C     gamma function  ratios  and their  inverse.   ACM  Trans.  Math.
84 C     Softw. 12 (1986), 377-393.
85 C
86 C     Computation of other parameters involve a seach for a value that
87 C     produces  the desired  value  of P.   The search relies  on  the
88 C     monotinicity of P with the other parameter.
89 C
90 C
91 C                              Note
92 C
93 C
94 C
95 C     The gamma density is proportional to
96 C       T**(SHAPE - 1) * EXP(- SCALE * T)
97 C
98 C                              History
99 C
100 C     Routine modified by Scilab group 1998, because of an undefined
101 C     variable. See below comments "CSS".
103 C**********************************************************************
104 C     .. Parameters ..
105       DOUBLE PRECISION tol
106       PARAMETER (tol=1.0D-13)
107       DOUBLE PRECISION atol
108       PARAMETER (atol=1.0D-50)
109       DOUBLE PRECISION zero,inf
110       PARAMETER (zero=1.0D-300,inf=1.0D300)
111 C     ..
112 C     .. Scalar Arguments ..
113       DOUBLE PRECISION bound,p,q,scale,shape,x
114       INTEGER status,which
115 C     ..
116 C     .. Local Scalars ..
117       DOUBLE PRECISION xx
118       DOUBLE PRECISION fx,xscale,cum,ccum,pq,porq
119       INTEGER ierr
120       LOGICAL qhi,qleft,qporq
121 C     ..
122 C     .. External Functions ..
123       DOUBLE PRECISION spmpar
124       EXTERNAL spmpar
125 C     ..
126 C     .. External Subroutines ..
127       EXTERNAL gaminv,dinvr,dstinv,cumgam
128 C     ..
129 C     .. Executable Statements ..
130 C
131 C     Check arguments
132 C
133       IF (.NOT. ((which.LT.1).OR. (which.GT.4))) GO TO 30
134       IF (.NOT. (which.LT.1)) GO TO 10
135       bound = 1.0D0
136       GO TO 20
138    10 bound = 4.0D0
139    20 status = -1
140       RETURN
142    30 IF (which.EQ.1) GO TO 70
143 C
144 C     P
145 C
146       IF (.NOT. ((p.LT.0.0D0).OR. (p.GT.1.0D0))) GO TO 60
147       IF (.NOT. (p.LT.0.0D0)) GO TO 40
148       bound = 0.0D0
149       GO TO 50
151    40 bound = 1.0d0
152    50 status = -2
153       RETURN
155    60 CONTINUE
156    70 IF (which.EQ.1) GO TO 110
157 C
158 C     Q
159 C
160       IF (.NOT. ((q.LE.0.0D0).OR. (q.GT.1.0D0))) GO TO 100
161       IF (.NOT. (q.LE.0.0D0)) GO TO 80
162       bound = 0.0D0
163       GO TO 90
165    80 bound = 1.0D0
166    90 status = -3
167       RETURN
169   100 CONTINUE
170   110 IF (which.EQ.2) GO TO 130
171 C
172 C     X
173 C
174       IF (.NOT. (x.LT.0.0D0)) GO TO 120
175       bound = 0.0D0
176       status = -4
177       RETURN
179   120 CONTINUE
180   130 IF (which.EQ.3) GO TO 150
181 C
182 C     SHAPE
183 C
184       IF (.NOT. (shape.LE.0.0D0)) GO TO 140
185       bound = 0.0D0
186       status = -5
187       RETURN
189   140 CONTINUE
190   150 IF (which.EQ.4) GO TO 170
191 C
192 C     SCALE
193 C
194       IF (.NOT. (scale.LE.0.0D0)) GO TO 160
195       bound = 0.0D0
196       status = -6
197       RETURN
199   160 CONTINUE
200   170 IF (which.EQ.1) GO TO 210
201 C
202 C     P + Q
203 C
204       pq = p + q
205       IF (.NOT. (abs(((pq)-0.5D0)-0.5D0).GT.
206      +    (3.0D0*spmpar(1)))) GO TO 200
207       IF (.NOT. (pq.LT.0.0D0)) GO TO 180
208       bound = 0.0D0
209       GO TO 190
211   180 bound = 1.0D0
212   190 status = 3
213       RETURN
215   200 CONTINUE
216   210 IF (which.EQ.1) GO TO 240
217 C
218 C     Select the minimum of P or Q
219 C
220       qporq = p .LE. q
221       IF (.NOT. (qporq)) GO TO 220
222       porq = p
223       GO TO 230
225   220 porq = q
226   230 CONTINUE
227 C
229 C
230   240 IF ((1).EQ. (which)) THEN
231 C
232 C     Calculating P
233 C
234           status = 0
235           xscale = x*scale
236           CALL cumgam(xscale,shape,p,q)
237 CSS   Next line changed by Scilab group. porq undefined here
238 CSS       IF (porq.GT.1.5D0) status = 10
239           IF (p.GT.1.5D0) status = 10
241       ELSE IF ((2).EQ. (which)) THEN
242 C
243 C     Computing X
244 C
245           CALL gaminv(shape,xx,-1.0D0,p,q,ierr)
246           IF (ierr.LT.0.0D0) THEN
247               status = 10
248               RETURN
250           ELSE
251               x = xx/scale
252               status = 0
253           END IF
255       ELSE IF ((3).EQ. (which)) THEN
256 C
257 C     Computing SHAPE
258 C
259           shape = 5.0D0
260           xscale = x*scale
261           CALL dstinv(zero,inf,0.5D0,0.5D0,5.0D0,atol,tol)
262           status = 0
263           CALL dinvr(status,shape,fx,qleft,qhi)
264   250     IF (.NOT. (status.EQ.1)) GO TO 290
265           CALL cumgam(xscale,shape,cum,ccum)
266           IF (.NOT. (qporq)) GO TO 260
267           fx = cum - p
268           GO TO 270
270   260     fx = ccum - q
271   270     IF (.NOT. ((qporq.AND. (cum.GT.1.5D0)).OR.
272      +        ((.NOT.qporq).AND. (ccum.GT.1.5D0)))) GO TO 280
273           status = 10
274           RETURN
276   280     CALL dinvr(status,shape,fx,qleft,qhi)
277           GO TO 250
279   290     IF (.NOT. (status.EQ.-1)) GO TO 320
280           IF (.NOT. (qleft)) GO TO 300
281           status = 1
282           bound = zero
283           GO TO 310
285   300     status = 2
286           bound = inf
287   310     CONTINUE
288   320     CONTINUE
290       ELSE IF ((4).EQ. (which)) THEN
291 C
292 C     Computing SCALE
293 C
294           CALL gaminv(shape,xx,-1.0D0,p,q,ierr)
295           IF (ierr.LT.0.0D0) THEN
296               status = 10
297               RETURN
299           ELSE
300               scale = xx/x
301               status = 0
302           END IF
304       END IF
306       RETURN
308       END