Statistics: fix for cdfnor 64bits
[scilab.git] / scilab / modules / statistics / src / dcdflib / cdfnor.f
1       SUBROUTINE cdfnor(which,p,q,x,mean,sd,status,bound)
2 C**********************************************************************
3 C
4 C      SUBROUTINE CDFNOR( WHICH, P, Q, X, MEAN, SD, STATUS, BOUND )
5 C               Cumulative Distribution Function
6 C               NORmal distribution
7 C
8 C
9 C                              Function
10 C
11 C
12 C     Calculates any one parameter of the normal
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  parameter
20 C     values is to be calculated using values  of the others.
21 C     Legal range: 1..4
22 C               iwhich = 1 : Calculate P and Q from X,MEAN and SD
23 C               iwhich = 2 : Calculate X from P,Q,MEAN and SD
24 C               iwhich = 3 : Calculate MEAN from P,Q,X and SD
25 C               iwhich = 4 : Calculate SD from P,Q,X and MEAN
26 C                    INTEGER WHICH
27 C
28 C     P <--> The integral from -infinity to X of the normal 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     X < --> Upper limit of integration of the normal-density.
38 C             Input range: ( -infinity, +infinity)
39 C                    DOUBLE PRECISION X
40 C
41 C     MEAN <--> The mean of the normal density.
42 C               Input range: (-infinity, +infinity)
43 C                    DOUBLE PRECISION MEAN
44 C
45 C     SD <--> Standard Deviation of the normal density.
46 C             Input range: (0, +infinity).
47 C                    DOUBLE PRECISION SD
48 C
49 C     STATUS <-- 0 if calculation completed correctly
50 C               -I if input parameter number I is out of range
51 C                1 if answer appears to be lower than lowest
52 C                  search bound
53 C                2 if answer appears to be higher than greatest
54 C                  search bound
55 C                3 if P + Q .ne. 1
56 C                    INTEGER STATUS
57 C
58 C     BOUND <-- Undefined if STATUS is 0
59 C
60 C               Bound exceeded by parameter number I if STATUS
61 C               is negative.
62 C
63 C               Lower search bound if STATUS is 1.
64 C
65 C               Upper search bound if STATUS is 2.
66 C
67 C
68 C                              Method
69 C
70 C
71 C
72 C
73 C     A slightly modified version of ANORM from
74 C
75 C     Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN
76 C     Package of Special Function Routines and Test Drivers"
77 C     acm Transactions on Mathematical Software. 19, 22-32.
78 C
79 C     is used to calulate the  cumulative standard normal distribution.
80 C
81 C     The rational functions from pages  90-95  of Kennedy and Gentle,
82 C     Statistical  Computing,  Marcel  Dekker, NY,  1980 are  used  as
83 C     starting values to Newton's Iterations which compute the inverse
84 C     standard normal.  Therefore no  searches  are necessary for  any
85 C     parameter.
86 C
87 C     For X < -15, the asymptotic expansion for the normal is used  as
88 C     the starting value in finding the inverse standard normal.
89 C     This is formula 26.2.12 of Abramowitz and Stegun.
90 C
91 C
92 C                              Note
93 C
94 C
95 C      The normal density is proportional to
96 C      exp( - 0.5 * (( X - MEAN)/SD)**2)
97 C
98 C
99 C**********************************************************************
100 C     .. Parameters ..
101       DOUBLE PRECISION inf
102       PARAMETER (inf=1.0D300)
103 C     ..
104 C     .. Scalar Arguments ..
105       DOUBLE PRECISION bound,mean,p,q,sd,x
106       INTEGER status,which
107 C     ..
108 C     .. Local Scalars ..
109       DOUBLE PRECISION z,pq
110 C     ..
111 C     .. External Functions ..
112       INTEGER vfinite
113       DOUBLE PRECISION dinvnr,spmpar
114       EXTERNAL dinvnr,spmpar
115 C     ..
116 C     .. External Subroutines ..
117       EXTERNAL cumnor
118 C     ..
119 C     .. Executable Statements ..
120 C
121 C     Check arguments
122 C
123       status = 0
124       IF (.NOT. ((which.LT.1).OR. (which.GT.4))) GO TO 30
125       IF (.NOT. (which.LT.1)) GO TO 10
126       bound = 1.0D0
127       GO TO 20
128
129    10 bound = 4.0D0
130    20 status = -1
131       RETURN
132
133    30 IF (which.EQ.1) GO TO 70
134 C
135 C     P
136 C
137       IF (ISANAN(p).EQ.1) THEN
138          CALL RETURNANANFORTRAN(x)
139          CALL RETURNANANFORTRAN(sd)
140          CALL RETURNANANFORTRAN(mean)
141          RETURN
142       ENDIF
143       IF (.NOT. ((p.LE.0.0D0).OR. (p.GT.1.0D0))) GO TO 60
144       IF (.NOT. (p.LE.0.0D0)) GO TO 40
145       bound = 0.0D0
146       GO TO 50
147
148    40 bound = 1.0D0
149    50 status = -2
150       RETURN
151
152    60 CONTINUE
153    70 IF (which.EQ.1) GO TO 110
154 C
155 C     Q
156 C
157       IF (ISANAN(q).EQ.1) THEN
158          CALL RETURNANANFORTRAN(x)
159          CALL RETURNANANFORTRAN(sd)
160          CALL RETURNANANFORTRAN(mean)
161          RETURN
162       ENDIF
163       IF (.NOT. ((q.LE.0.0D0).OR. (q.GT.1.0D0))) GO TO 100
164       IF (.NOT. (q.LE.0.0D0)) GO TO 80
165       bound = 0.0D0
166       GO TO 90
167
168    80 bound = 1.0D0
169    90 status = -3
170       RETURN
171
172   100 CONTINUE
173   110 IF (which.EQ.1) GO TO 150
174 C
175 C     P + Q
176 C
177       pq = p + q
178       IF (.NOT. (abs(((pq)-0.5D0)-0.5D0).GT.
179      +    (3.0D0*spmpar(1)))) GO TO 140
180       IF (.NOT. (pq.LT.0.0D0)) GO TO 120
181       bound = 0.0D0
182       GO TO 130
183
184   120 bound = 1.0D0
185   130 status = 3
186       RETURN
187
188   140 CONTINUE
189   150 IF (which.EQ.4) GO TO 170
190 C
191 C     X
192 C
193       IF (ISANAN(x).EQ.1) THEN
194          CALL RETURNANANFORTRAN(p)
195          CALL RETURNANANFORTRAN(q)
196          CALL RETURNANANFORTRAN(mean)
197          CALL RETURNANANFORTRAN(sd)
198          RETURN
199       ENDIF
200       IF (vfinite(1,x).EQ.0) then
201          IF (which.EQ.1) then
202             IF (x.GT.0) then
203                p = 1
204                q = 0
205                RETURN
206             ELSE
207                p = 0
208                q = 1
209                RETURN
210             ENDIF
211          ELSE
212             x = SIGN(inf,x)
213          ENDIF
214       ENDIF
215 C
216 C     MEAN
217 C
218       IF (ISANAN(mean).EQ.1) THEN
219          CALL RETURNANANFORTRAN(p)
220          CALL RETURNANANFORTRAN(q)
221          CALL RETURNANANFORTRAN(x)
222          CALL RETURNANANFORTRAN(sd)
223          RETURN
224       ENDIF
225       IF (vfinite(1,mean).EQ.0) mean = SIGN(inf,mean)
226 C
227 C     SD
228 C
229       IF (ISANAN(sd).EQ.1) THEN
230          CALL RETURNANANFORTRAN(p)
231          CALL RETURNANANFORTRAN(q)
232          CALL RETURNANANFORTRAN(x)
233          CALL RETURNANANFORTRAN(mean)
234          RETURN
235       ENDIF
236       IF (vfinite(1,sd).EQ.0) sd = SIGN(inf,sd)
237       IF (.NOT. (sd.LE.0.0D0)) GO TO 160
238       bound = 0.0D0
239       status = -6
240       RETURN
241
242   160 CONTINUE
243 C
244 C     Calculate ANSWERS
245 C
246   170 IF ((1).EQ. (which)) THEN
247 C
248 C     Computing P
249 C
250          z = (x-mean)/sd
251          CALL cumnor(z,p,q)
252
253       ELSE IF ((2).EQ. (which)) THEN
254 C
255 C     Computing X
256 C
257           z = dinvnr(p,q)
258           x = sd*z + mean
259
260       ELSE IF ((3).EQ. (which)) THEN
261 C
262 C     Computing the MEAN
263 C
264           z = dinvnr(p,q)
265           mean = x - sd*z
266
267       ELSE IF ((4).EQ. (which)) THEN
268 C
269 C     Computing SD
270 C
271           z = dinvnr(p,q)
272           sd = (x-mean)/z
273       END IF
274
275       RETURN
276
277       END