d0a73895643b941bf97e02640fff7850be798afb
[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 C     ..
102 C     .. Scalar Arguments ..
103       DOUBLE PRECISION bound,mean,p,q,sd,x
104       INTEGER status,which
105 C     ..
106 C     .. Local Scalars ..
107       DOUBLE PRECISION z,pq
108 C     ..
109 C     .. External Functions ..
110       INTEGER vfinite
111       DOUBLE PRECISION dinvnr,spmpar
112       EXTERNAL dinvnr,spmpar
113 C     ..
114 C     .. External Subroutines ..
115       EXTERNAL cumnor
116 C     ..
117 C     .. Executable Statements ..
118 C
119 C     Check arguments
120 C
121       status = 0
122       IF (.NOT. ((which.LT.1).OR. (which.GT.4))) GO TO 30
123       IF (.NOT. (which.LT.1)) GO TO 10
124       bound = 1.0D0
125       GO TO 20
126
127    10 bound = 4.0D0
128    20 status = -1
129       RETURN
130
131    30 IF (which.EQ.1) GO TO 70
132 C
133 C     P
134 C
135       IF (ISANAN(p).EQ.1) THEN
136          CALL RETURNANANFORTRAN(x)
137          CALL RETURNANANFORTRAN(sd)
138          CALL RETURNANANFORTRAN(mean)
139          RETURN
140       ENDIF
141       IF (.NOT. ((p.LE.0.0D0).OR. (p.GT.1.0D0))) GO TO 60
142       IF (.NOT. (p.LE.0.0D0)) GO TO 40
143       bound = 0.0D0
144       GO TO 50
145
146    40 bound = 1.0D0
147    50 status = -2
148       RETURN
149
150    60 CONTINUE
151    70 IF (which.EQ.1) GO TO 110
152 C
153 C     Q
154 C
155       IF (ISANAN(q).EQ.1) THEN
156          CALL RETURNANANFORTRAN(x)
157          CALL RETURNANANFORTRAN(sd)
158          CALL RETURNANANFORTRAN(mean)
159          RETURN
160       ENDIF
161       IF (.NOT. ((q.LE.0.0D0).OR. (q.GT.1.0D0))) GO TO 100
162       IF (.NOT. (q.LE.0.0D0)) GO TO 80
163       bound = 0.0D0
164       GO TO 90
165
166    80 bound = 1.0D0
167    90 status = -3
168       RETURN
169
170   100 CONTINUE
171   110 IF (which.EQ.1) GO TO 150
172 C
173 C     P + Q
174 C
175       pq = p + q
176       IF (.NOT. (abs(((pq)-0.5D0)-0.5D0).GT.
177      +    (3.0D0*spmpar(1)))) GO TO 140
178       IF (.NOT. (pq.LT.0.0D0)) GO TO 120
179       bound = 0.0D0
180       GO TO 130
181
182   120 bound = 1.0D0
183   130 status = 3
184       RETURN
185
186   140 CONTINUE
187   150 IF (which.EQ.4) GO TO 170
188 C
189 C     X
190 C
191       IF (ISANAN(x).EQ.1) THEN
192          CALL RETURNANANFORTRAN(p)
193          CALL RETURNANANFORTRAN(q)
194          CALL RETURNANANFORTRAN(mean)
195          CALL RETURNANANFORTRAN(sd)
196          RETURN
197       ENDIF
198       IF (vfinite(1,x).EQ.0) then
199          IF (which.EQ.1) then
200             IF (x.GT.0) then
201                p = 1
202                q = 0
203                RETURN
204             ELSE
205                p = 0
206                q = 1
207                RETURN
208             ENDIF
209          ELSE
210             x = SIGN(1D308,x)
211          ENDIF
212       ENDIF
213 C
214 C     MEAN
215 C
216       IF (ISANAN(mean).EQ.1) THEN
217          CALL RETURNANANFORTRAN(p)
218          CALL RETURNANANFORTRAN(q)
219          CALL RETURNANANFORTRAN(x)
220          CALL RETURNANANFORTRAN(sd)
221          RETURN
222       ENDIF
223       IF (vfinite(1,mean).EQ.0) mean = SIGN(1D308,mean)
224 C
225 C     SD
226 C
227       IF (ISANAN(sd).EQ.1) THEN
228          CALL RETURNANANFORTRAN(p)
229          CALL RETURNANANFORTRAN(q)
230          CALL RETURNANANFORTRAN(x)
231          CALL RETURNANANFORTRAN(mean)
232          RETURN
233       ENDIF
234       IF (vfinite(1,sd).EQ.0) sd = SIGN(1D308,sd)
235       IF (.NOT. (sd.LE.0.0D0)) GO TO 160
236       bound = 0.0D0
237       status = -6
238       RETURN
239
240   160 CONTINUE
241 C
242 C     Calculate ANSWERS
243 C
244   170 IF ((1).EQ. (which)) THEN
245 C
246 C     Computing P
247 C
248          z = (x-mean)/sd
249          CALL cumnor(z,p,q)
250
251       ELSE IF ((2).EQ. (which)) THEN
252 C
253 C     Computing X
254 C
255           z = dinvnr(p,q)
256           x = sd*z + mean
257
258       ELSE IF ((3).EQ. (which)) THEN
259 C
260 C     Computing the MEAN
261 C
262           z = dinvnr(p,q)
263           mean = x - sd*z
264
265       ELSE IF ((4).EQ. (which)) THEN
266 C
267 C     Computing SD
268 C
269           z = dinvnr(p,q)
270           sd = (x-mean)/z
271       END IF
272
273       RETURN
274
275       END