* Bug #7296 fixed - Enabling %nan, %inf and -%inf for the cdf* functions
[scilab.git] / scilab / modules / statistics / src / dcdflib / cdfgam.f
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".
102
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       INTEGER vfinite
124       DOUBLE PRECISION spmpar
125       EXTERNAL spmpar
126 C     ..
127 C     .. External Subroutines ..
128       EXTERNAL gaminv,dinvr,dstinv,cumgam
129 C     ..
130 C     .. Executable Statements ..
131 C
132 C     Check arguments
133 C
134       IF (.NOT. ((which.LT.1).OR. (which.GT.4))) GO TO 30
135       IF (.NOT. (which.LT.1)) GO TO 10
136       bound = 1.0D0
137       GO TO 20
138
139    10 bound = 4.0D0
140    20 status = -1
141       RETURN
142
143    30 IF (which.EQ.1) GO TO 70
144 C
145 C     P
146 C
147       IF (ISANAN(p).EQ.1) THEN
148          CALL RETURNANANFORTRAN(shape)
149          CALL RETURNANANFORTRAN(x)
150          CALL RETURNANANFORTRAN(scale)
151          RETURN
152       ENDIF
153       IF (.NOT. ((p.LT.0.0D0).OR. (p.GT.1.0D0))) GO TO 60
154       IF (.NOT. (p.LT.0.0D0)) GO TO 40
155       bound = 0.0D0
156       GO TO 50
157
158    40 bound = 1.0d0
159    50 status = -2
160       RETURN
161
162    60 CONTINUE
163    70 IF (which.EQ.1) GO TO 110
164 C
165 C     Q
166 C
167       IF (ISANAN(q).EQ.1) THEN
168          CALL RETURNANANFORTRAN(shape)
169          CALL RETURNANANFORTRAN(x)
170          CALL RETURNANANFORTRAN(scale)
171          RETURN
172       ENDIF
173       IF (.NOT. ((q.LE.0.0D0).OR. (q.GT.1.0D0))) GO TO 100
174       IF (.NOT. (q.LE.0.0D0)) GO TO 80
175       bound = 0.0D0
176       GO TO 90
177
178    80 bound = 1.0D0
179    90 status = -3
180       RETURN
181
182   100 CONTINUE
183   110 IF (which.EQ.2) GO TO 130
184 C
185 C     X
186 C
187       IF (ISANAN(x).EQ.1) THEN
188          CALL RETURNANANFORTRAN(p)
189          CALL RETURNANANFORTRAN(q)
190          CALL RETURNANANFORTRAN(shape)
191          CALL RETURNANANFORTRAN(scale)
192          RETURN
193       ENDIF
194       IF (vfinite(1,x).EQ.0) then
195          IF (which.EQ.1) then
196             IF (x.GT.0) then
197                p = 1
198                q = 0
199                RETURN
200             ENDIF
201          ELSE
202             x = SIGN(1D300,x)
203          ENDIF
204       ENDIF
205       IF (.NOT. (x.LT.0.0D0)) GO TO 120
206       bound = 0.0D0
207       status = -4
208       RETURN
209
210   120 CONTINUE
211   130 IF (which.EQ.3) GO TO 150
212 C
213 C     SHAPE
214 C
215       IF (ISANAN(shape).EQ.1) THEN
216          CALL RETURNANANFORTRAN(p)
217          CALL RETURNANANFORTRAN(q)
218          CALL RETURNANANFORTRAN(x)
219          CALL RETURNANANFORTRAN(scale)
220          RETURN
221       ENDIF
222       IF (vfinite(1,shape).EQ.0) shape = SIGN(inf,shape)
223       IF (.NOT. (shape.LE.0.0D0)) GO TO 140
224       bound = 0.0D0
225       status = -5
226       RETURN
227
228   140 CONTINUE
229   150 IF (which.EQ.4) GO TO 170
230 C
231 C     SCALE
232 C
233       IF (ISANAN(scale).EQ.1) THEN
234          CALL RETURNANANFORTRAN(p)
235          CALL RETURNANANFORTRAN(q)
236          CALL RETURNANANFORTRAN(shape)
237          CALL RETURNANANFORTRAN(x)
238          RETURN
239       ENDIF
240       IF (vfinite(1,scale).EQ.0) scale = SIGN(inf,scale)
241       IF (.NOT. (scale.LE.0.0D0)) GO TO 160
242       bound = 0.0D0
243       status = -6
244       RETURN
245
246   160 CONTINUE
247   170 IF (which.EQ.1) GO TO 210
248 C
249 C     P + Q
250 C
251       pq = p + q
252       IF (.NOT. (abs(((pq)-0.5D0)-0.5D0).GT.
253      +    (3.0D0*spmpar(1)))) GO TO 200
254       IF (.NOT. (pq.LT.0.0D0)) GO TO 180
255       bound = 0.0D0
256       GO TO 190
257
258   180 bound = 1.0D0
259   190 status = 3
260       RETURN
261
262   200 CONTINUE
263   210 IF (which.EQ.1) GO TO 240
264 C
265 C     Select the minimum of P or Q
266 C
267       qporq = p .LE. q
268       IF (.NOT. (qporq)) GO TO 220
269       porq = p
270       GO TO 230
271
272   220 porq = q
273   230 CONTINUE
274 C
275 C     Calculate ANSWERS
276 C
277   240 IF ((1).EQ. (which)) THEN
278 C
279 C     Calculating P
280 C
281           status = 0
282           xscale = x*scale
283           CALL cumgam(xscale,shape,p,q)
284 CSS   Next line changed by Scilab group. porq undefined here
285 CSS       IF (porq.GT.1.5D0) status = 10
286           IF (p.GT.1.5D0) status = 10
287
288       ELSE IF ((2).EQ. (which)) THEN
289 C
290 C     Computing X
291 C
292           CALL gaminv(shape,xx,-1.0D0,p,q,ierr)
293           IF (ierr.LT.0.0D0) THEN
294               status = 10
295               RETURN
296
297           ELSE
298               x = xx/scale
299               status = 0
300           END IF
301
302       ELSE IF ((3).EQ. (which)) THEN
303 C
304 C     Computing SHAPE
305 C
306           shape = 5.0D0
307           xscale = x*scale
308           CALL dstinv(zero,inf,0.5D0,0.5D0,5.0D0,atol,tol)
309           status = 0
310           CALL dinvr(status,shape,fx,qleft,qhi)
311   250     IF (.NOT. (status.EQ.1)) GO TO 290
312           CALL cumgam(xscale,shape,cum,ccum)
313           IF (.NOT. (qporq)) GO TO 260
314           fx = cum - p
315           GO TO 270
316
317   260     fx = ccum - q
318   270     IF (.NOT. ((qporq.AND. (cum.GT.1.5D0)).OR.
319      +        ((.NOT.qporq).AND. (ccum.GT.1.5D0)))) GO TO 280
320           status = 10
321           RETURN
322
323   280     CALL dinvr(status,shape,fx,qleft,qhi)
324           GO TO 250
325
326   290     IF (.NOT. (status.EQ.-1)) GO TO 320
327           IF (.NOT. (qleft)) GO TO 300
328           status = 1
329           bound = zero
330           GO TO 310
331
332   300     status = 2
333           bound = inf
334   310     CONTINUE
335   320     CONTINUE
336
337       ELSE IF ((4).EQ. (which)) THEN
338 C
339 C     Computing SCALE
340 C
341           CALL gaminv(shape,xx,-1.0D0,p,q,ierr)
342           IF (ierr.LT.0.0D0) THEN
343               status = 10
344               RETURN
345
346           ELSE
347               scale = xx/x
348               status = 0
349           END IF
350
351       END IF
352
353       RETURN
354
355       END