← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to scipy/special/cdflib/cdfnor.f

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
scipy/special/cdflib/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     .. Scalar Arguments ..
 
101
      DOUBLE PRECISION bound,mean,p,q,sd,x
 
102
      INTEGER status,which
 
103
C     ..
 
104
C     .. Local Scalars ..
 
105
      DOUBLE PRECISION pq,z
 
106
C     ..
 
107
C     .. External Functions ..
 
108
      DOUBLE PRECISION dinvnr,spmpar
 
109
      EXTERNAL dinvnr,spmpar
 
110
C     ..
 
111
C     .. External Subroutines ..
 
112
      EXTERNAL cumnor
 
113
C     ..
 
114
C     .. Intrinsic Functions ..
 
115
      INTRINSIC abs
 
116
C     ..
 
117
      status = 0
 
118
      IF (.NOT. ((which.LT.1).OR. (which.GT.4))) GO TO 30
 
119
      IF (.NOT. (which.LT.1)) GO TO 10
 
120
      bound = 1.0D0
 
121
      GO TO 20
 
122
 
 
123
   10 bound = 4.0D0
 
124
   20 status = -1
 
125
      RETURN
 
126
 
 
127
   30 IF (which.EQ.1) GO TO 70
 
128
      IF (.NOT. ((p.LE.0.0D0).OR. (p.GT.1.0D0))) GO TO 60
 
129
      IF (.NOT. (p.LE.0.0D0)) GO TO 40
 
130
      bound = 0.0D0
 
131
      GO TO 50
 
132
 
 
133
   40 bound = 1.0D0
 
134
   50 status = -2
 
135
      RETURN
 
136
 
 
137
   60 CONTINUE
 
138
   70 IF (which.EQ.1) GO TO 110
 
139
      IF (.NOT. ((q.LE.0.0D0).OR. (q.GT.1.0D0))) GO TO 100
 
140
      IF (.NOT. (q.LE.0.0D0)) GO TO 80
 
141
      bound = 0.0D0
 
142
      GO TO 90
 
143
 
 
144
   80 bound = 1.0D0
 
145
   90 status = -3
 
146
      RETURN
 
147
 
 
148
  100 CONTINUE
 
149
  110 IF (which.EQ.1) GO TO 150
 
150
      pq = p + q
 
151
      IF (.NOT. (abs(((pq)-0.5D0)-0.5D0).GT.
 
152
     +    (3.0D0*spmpar(1)))) GO TO 140
 
153
      IF (.NOT. (pq.LT.0.0D0)) GO TO 120
 
154
      bound = 0.0D0
 
155
      GO TO 130
 
156
 
 
157
  120 bound = 1.0D0
 
158
  130 status = 3
 
159
      RETURN
 
160
 
 
161
  140 CONTINUE
 
162
  150 IF (which.EQ.4) GO TO 170
 
163
      IF (.NOT. (sd.LE.0.0D0)) GO TO 160
 
164
      bound = 0.0D0
 
165
      status = -6
 
166
      RETURN
 
167
 
 
168
  160 CONTINUE
 
169
  170 IF ((1).EQ. (which)) THEN
 
170
          z = (x-mean)/sd
 
171
          CALL cumnor(z,p,q)
 
172
 
 
173
      ELSE IF ((2).EQ. (which)) THEN
 
174
          z = dinvnr(p,q)
 
175
          x = sd*z + mean
 
176
 
 
177
      ELSE IF ((3).EQ. (which)) THEN
 
178
          z = dinvnr(p,q)
 
179
          mean = x - sd*z
 
180
 
 
181
      ELSE IF ((4).EQ. (which)) THEN
 
182
          z = dinvnr(p,q)
 
183
          sd = (x-mean)/z
 
184
      END IF
 
185
 
 
186
      RETURN
 
187
 
 
188
      END