← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to routines/dcd/cdfpoi.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/dcd/cdfpoi.f
 
1
      SUBROUTINE cdfpoi(which,p,q,s,xlam,status,bound)
 
2
C**********************************************************************
 
3
C
 
4
C      SUBROUTINE CDFPOI( WHICH, P, Q, S, XLAM, STATUS, BOUND )
 
5
C               Cumulative Distribution Function
 
6
C               POIsson distribution
 
7
C
 
8
C
 
9
C                              Function
 
10
C
 
11
C
 
12
C     Calculates any one parameter of the Poisson
 
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  argument
 
20
C               value is to be calculated from the others.
 
21
C               Legal range: 1..3
 
22
C               iwhich = 1 : Calculate P and Q from S and XLAM
 
23
C               iwhich = 2 : Calculate A from P,Q and XLAM
 
24
C               iwhich = 3 : Calculate XLAM from P,Q and S
 
25
C                    INTEGER WHICH
 
26
C
 
27
C        P <--> The cumulation from 0 to S of the poisson density.
 
28
C               Input range: [0,1].
 
29
C                    DOUBLE PRECISION P
 
30
C
 
31
C        Q <--> 1-P.
 
32
C               Input range: (0, 1].
 
33
C               P + Q = 1.0.
 
34
C                    DOUBLE PRECISION Q
 
35
C
 
36
C        S <--> Upper limit of cumulation of the Poisson.
 
37
C               Input range: [0, +infinity).
 
38
C               Search range: [0,1E300]
 
39
C                    DOUBLE PRECISION S
 
40
C
 
41
C     XLAM <--> Mean of the Poisson distribution.
 
42
C               Input range: [0, +infinity).
 
43
C               Search range: [0,1E300]
 
44
C                    DOUBLE PRECISION XLAM
 
45
C
 
46
C     STATUS <-- 0 if calculation completed correctly
 
47
C               -I if input parameter number I is out of range
 
48
C                1 if answer appears to be lower than lowest
 
49
C                  search bound
 
50
C                2 if answer appears to be higher than greatest
 
51
C                  search bound
 
52
C                3 if P + Q .ne. 1
 
53
C                    INTEGER STATUS
 
54
C
 
55
C     BOUND <-- Undefined if STATUS is 0
 
56
C
 
57
C               Bound exceeded by parameter number I if STATUS
 
58
C               is negative.
 
59
C
 
60
C               Lower search bound if STATUS is 1.
 
61
C
 
62
C               Upper search bound if STATUS is 2.
 
63
C
 
64
C
 
65
C                              Method
 
66
C
 
67
C
 
68
C     Formula   26.4.21  of   Abramowitz  and   Stegun,   Handbook  of
 
69
C     Mathematical Functions (1966) is used  to reduce the computation
 
70
C     of  the cumulative distribution function to that  of computing a
 
71
C     chi-square, hence an incomplete gamma function.
 
72
C
 
73
C     Cumulative  distribution function  (P) is  calculated  directly.
 
74
C     Computation of other parameters involve a seach for a value that
 
75
C     produces  the desired value of  P.   The  search relies  on  the
 
76
C     monotinicity of P with the other parameter.
 
77
C
 
78
C
 
79
C**********************************************************************
 
80
C     .. Parameters ..
 
81
      DOUBLE PRECISION tol
 
82
      PARAMETER (tol=1.0D-8)
 
83
      DOUBLE PRECISION atol
 
84
      PARAMETER (atol=1.0D-50)
 
85
      DOUBLE PRECISION inf
 
86
      PARAMETER (inf=1.0D300)
 
87
C     ..
 
88
C     .. Scalar Arguments ..
 
89
      DOUBLE PRECISION bound,p,q,s,xlam
 
90
      INTEGER status,which
 
91
C     ..
 
92
C     .. Local Scalars ..
 
93
      DOUBLE PRECISION fx,cum,ccum,pq
 
94
      LOGICAL qhi,qleft,qporq
 
95
C     ..
 
96
C     .. External Functions ..
 
97
      DOUBLE PRECISION spmpar
 
98
      EXTERNAL spmpar
 
99
C     ..
 
100
C     .. External Subroutines ..
 
101
      EXTERNAL dinvr,dstinv,cumpoi
 
102
C     ..
 
103
C     .. Executable Statements ..
 
104
C
 
105
C     Check arguments
 
106
C
 
107
      IF (.NOT. ((which.LT.1).OR. (which.GT.3))) GO TO 30
 
108
      IF (.NOT. (which.LT.1)) GO TO 10
 
109
      bound = 1.0D0
 
110
      GO TO 20
 
111
 
 
112
   10 bound = 3.0D0
 
113
   20 status = -1
 
114
      RETURN
 
115
 
 
116
   30 IF (which.EQ.1) GO TO 70
 
117
C
 
118
C     P
 
119
C
 
120
      IF (.NOT. ((p.LT.0.0D0).OR. (p.GT.1.0D0))) GO TO 60
 
121
      IF (.NOT. (p.LT.0.0D0)) GO TO 40
 
122
      bound = 0.0D0
 
123
      GO TO 50
 
124
 
 
125
   40 bound = 1.0D0
 
126
   50 status = -2
 
127
      RETURN
 
128
 
 
129
   60 CONTINUE
 
130
   70 IF (which.EQ.1) GO TO 110
 
131
C
 
132
C     Q
 
133
C
 
134
      IF (.NOT. ((q.LE.0.0D0).OR. (q.GT.1.0D0))) GO TO 100
 
135
      IF (.NOT. (q.LE.0.0D0)) GO TO 80
 
136
      bound = 0.0D0
 
137
      GO TO 90
 
138
 
 
139
   80 bound = 1.0D0
 
140
   90 status = -3
 
141
      RETURN
 
142
 
 
143
  100 CONTINUE
 
144
  110 IF (which.EQ.2) GO TO 130
 
145
C
 
146
C     S
 
147
C
 
148
      IF (.NOT. (s.LT.0.0D0)) GO TO 120
 
149
      bound = 0.0D0
 
150
      status = -4
 
151
      RETURN
 
152
 
 
153
  120 CONTINUE
 
154
  130 IF (which.EQ.3) GO TO 150
 
155
C
 
156
C     XLAM
 
157
C
 
158
      IF (.NOT. (xlam.LT.0.0D0)) GO TO 140
 
159
      bound = 0.0D0
 
160
      status = -5
 
161
      RETURN
 
162
 
 
163
  140 CONTINUE
 
164
  150 IF (which.EQ.1) GO TO 190
 
165
C
 
166
C     P + Q
 
167
C
 
168
      pq = p + q
 
169
      IF (.NOT. (abs(((pq)-0.5D0)-0.5D0).GT.
 
170
     +    (3.0D0*spmpar(1)))) GO TO 180
 
171
      IF (.NOT. (pq.LT.0.0D0)) GO TO 160
 
172
      bound = 0.0D0
 
173
      GO TO 170
 
174
 
 
175
  160 bound = 1.0D0
 
176
  170 status = 3
 
177
      RETURN
 
178
 
 
179
  180 CONTINUE
 
180
  190 IF (.NOT. (which.EQ.1)) qporq = p .LE. q
 
181
C
 
182
C     Select the minimum of P or Q
 
183
C
 
184
C
 
185
C     Calculate ANSWERS
 
186
C
 
187
      IF ((1).EQ. (which)) THEN
 
188
C
 
189
C     Calculating P
 
190
C
 
191
          CALL cumpoi(s,xlam,p,q)
 
192
          status = 0
 
193
 
 
194
      ELSE IF ((2).EQ. (which)) THEN
 
195
C
 
196
C     Calculating S
 
197
C
 
198
          s = 5.0D0
 
199
          CALL dstinv(0.0D0,inf,0.5D0,0.5D0,5.0D0,atol,tol)
 
200
          status = 0
 
201
          CALL dinvr(status,s,fx,qleft,qhi)
 
202
  200     IF (.NOT. (status.EQ.1)) GO TO 230
 
203
          CALL cumpoi(s,xlam,cum,ccum)
 
204
          IF (.NOT. (qporq)) GO TO 210
 
205
          fx = cum - p
 
206
          GO TO 220
 
207
 
 
208
  210     fx = ccum - q
 
209
  220     CALL dinvr(status,s,fx,qleft,qhi)
 
210
          GO TO 200
 
211
 
 
212
  230     IF (.NOT. (status.EQ.-1)) GO TO 260
 
213
          IF (.NOT. (qleft)) GO TO 240
 
214
          status = 1
 
215
          bound = 0.0D0
 
216
          GO TO 250
 
217
 
 
218
  240     status = 2
 
219
          bound = inf
 
220
  250     CONTINUE
 
221
  260     CONTINUE
 
222
 
 
223
      ELSE IF ((3).EQ. (which)) THEN
 
224
C
 
225
C     Calculating XLAM
 
226
C
 
227
          xlam = 5.0D0
 
228
          CALL dstinv(0.0D0,inf,0.5D0,0.5D0,5.0D0,atol,tol)
 
229
          status = 0
 
230
          CALL dinvr(status,xlam,fx,qleft,qhi)
 
231
  270     IF (.NOT. (status.EQ.1)) GO TO 300
 
232
          CALL cumpoi(s,xlam,cum,ccum)
 
233
          IF (.NOT. (qporq)) GO TO 280
 
234
          fx = cum - p
 
235
          GO TO 290
 
236
 
 
237
  280     fx = ccum - q
 
238
  290     CALL dinvr(status,xlam,fx,qleft,qhi)
 
239
          GO TO 270
 
240
 
 
241
  300     IF (.NOT. (status.EQ.-1)) GO TO 330
 
242
          IF (.NOT. (qleft)) GO TO 310
 
243
          status = 1
 
244
          bound = 0.0D0
 
245
          GO TO 320
 
246
 
 
247
  310     status = 2
 
248
          bound = inf
 
249
  320     CONTINUE
 
250
  330 END IF
 
251
 
 
252
      RETURN
 
253
 
 
254
      END