1 /*
2 * Mathlib : A C Library of Special Functions
3 * Copyright (C) 1998 Ross Ihaka
4 * Copyright (C) 2000-2010 The R Core Team
5 * Copyright (C) 2003 The R Foundation
6 *
7 * This program is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program; if not, a copy is available at
19 * http://www.r-project.org/Licenses/
20 *
21 * SYNOPSIS
22 *
23 * #include <Rmath.h>
24 *
25 * double pnorm5(double x, double mu, double sigma, int lower_tail,int log_p);
26 * {pnorm (..) is synonymous and preferred inside R}
27 *
28 * void pnorm_both(double x, double *cum, double *ccum,
29 * int i_tail, int log_p);
30 *
31 * DESCRIPTION
32 *
33 * The main computation evaluates near-minimax approximations derived
34 * from those in "Rational Chebyshev approximations for the error
35 * function" by W. J. Cody, Math. Comp., 1969, 631-637. This
36 * transportable program uses rational functions that theoretically
37 * approximate the normal distribution function to at least 18
38 * significant decimal digits. The accuracy achieved depends on the
39 * arithmetic system, the compiler, the intrinsic functions, and
40 * proper selection of the machine-dependent constants.
41 *
42 * REFERENCE
43 *
44 * Cody, W. D. (1993).
45 * ALGORITHM 715: SPECFUN - A Portable FORTRAN Package of
46 * Special Function Routines and Test Drivers".
47 * ACM Transactions on Mathematical Software. 19, 22-32.
48 *
49 * EXTENSIONS
50 *
51 * The "_both" , lower, upper, and log_p variants were added by
52 * Martin Maechler, Jan.2000;
53 * as well as log1p() and similar improvements later on.
54 *
55 * James M. Rath contributed bug report PR#699 and patches correcting SIXTEN
56 * and if() clauses {with a bug: "|| instead of &&" -> PR #2883) more in line
57 * with the original Cody code.
58 */
59
60 #include "nmath.h"
61 #include "dpq.h"
pnorm5(double x,double mu,double sigma,int lower_tail,int log_p)62 double pnorm5(double x, double mu, double sigma, int lower_tail, int log_p)
63 {
64 double p, cp;
65
66 /* Note: The structure of these checks has been carefully thought through.
67 * For example, if x == mu and sigma == 0, we get the correct answer 1.
68 */
69 #ifdef IEEE_754
70 if(ISNAN(x) || ISNAN(mu) || ISNAN(sigma))
71 return x + mu + sigma;
72 #endif
73 if(!R_FINITE(x) && mu == x) return ML_NAN;/* x-mu is NaN */
74 if (sigma <= 0) {
75 if(sigma < 0) ML_ERR_return_NAN;
76 /* sigma = 0 : */
77 return (x < mu) ? R_DT_0 : R_DT_1;
78 }
79 p = (x - mu) / sigma;
80 if(!R_FINITE(p))
81 return (x < mu) ? R_DT_0 : R_DT_1;
82 x = p;
83
84 pnorm_both(x, &p, &cp, (lower_tail ? 0 : 1), log_p);
85
86 return(lower_tail ? p : cp);
87 }
88
89 #define SIXTEN 16 /* Cutoff allowing exact "*" and "/" */
90
pnorm_both(double x,double * cum,double * ccum,int i_tail,int log_p)91 void pnorm_both(double x, double *cum, double *ccum, int i_tail, int log_p)
92 {
93 /* i_tail in {0,1,2} means: "lower", "upper", or "both" :
94 if(lower) return *cum := P[X <= x]
95 if(upper) return *ccum := P[X > x] = 1 - P[X <= x]
96 */
97 const static double a[5] = {
98 2.2352520354606839287,
99 161.02823106855587881,
100 1067.6894854603709582,
101 18154.981253343561249,
102 0.065682337918207449113
103 };
104 const static double b[4] = {
105 47.20258190468824187,
106 976.09855173777669322,
107 10260.932208618978205,
108 45507.789335026729956
109 };
110 const static double c[9] = {
111 0.39894151208813466764,
112 8.8831497943883759412,
113 93.506656132177855979,
114 597.27027639480026226,
115 2494.5375852903726711,
116 6848.1904505362823326,
117 11602.651437647350124,
118 9842.7148383839780218,
119 1.0765576773720192317e-8
120 };
121 const static double d[8] = {
122 22.266688044328115691,
123 235.38790178262499861,
124 1519.377599407554805,
125 6485.558298266760755,
126 18615.571640885098091,
127 34900.952721145977266,
128 38912.003286093271411,
129 19685.429676859990727
130 };
131 const static double p[6] = {
132 0.21589853405795699,
133 0.1274011611602473639,
134 0.022235277870649807,
135 0.001421619193227893466,
136 2.9112874951168792e-5,
137 0.02307344176494017303
138 };
139 const static double q[5] = {
140 1.28426009614491121,
141 0.468238212480865118,
142 0.0659881378689285515,
143 0.00378239633202758244,
144 7.29751555083966205e-5
145 };
146
147 double xden, xnum, temp, del, eps, xsq, y;
148 #ifdef NO_DENORMS
149 double min = DBL_MIN;
150 #endif
151 int i, lower, upper;
152
153 #ifdef IEEE_754
154 if(ISNAN(x)) { *cum = *ccum = x; return; }
155 #endif
156
157 /* Consider changing these : */
158 eps = DBL_EPSILON * 0.5;
159
160 /* i_tail in {0,1,2} =^= {lower, upper, both} */
161 lower = i_tail != 1;
162 upper = i_tail != 0;
163
164 y = fabs(x);
165 if (y <= 0.67448975) { /* qnorm(3/4) = .6744.... -- earlier had 0.66291 */
166 if (y > eps) {
167 xsq = x * x;
168 xnum = a[4] * xsq;
169 xden = xsq;
170 for (i = 0; i < 3; ++i) {
171 xnum = (xnum + a[i]) * xsq;
172 xden = (xden + b[i]) * xsq;
173 }
174 } else xnum = xden = 0.0;
175
176 temp = x * (xnum + a[3]) / (xden + b[3]);
177 if(lower) *cum = 0.5 + temp;
178 if(upper) *ccum = 0.5 - temp;
179 if(log_p) {
180 if(lower) *cum = log(*cum);
181 if(upper) *ccum = log(*ccum);
182 }
183 }
184 else if (y <= M_SQRT_32) {
185
186 /* Evaluate pnorm for 0.674.. = qnorm(3/4) < |x| <= sqrt(32) ~= 5.657 */
187
188 xnum = c[8] * y;
189 xden = y;
190 for (i = 0; i < 7; ++i) {
191 xnum = (xnum + c[i]) * y;
192 xden = (xden + d[i]) * y;
193 }
194 temp = (xnum + c[7]) / (xden + d[7]);
195
196 #define do_del(X) \
197 xsq = trunc(X * SIXTEN) / SIXTEN; \
198 del = (X - xsq) * (X + xsq); \
199 if(log_p) { \
200 *cum = (-xsq * xsq * 0.5) + (-del * 0.5) + log(temp); \
201 if((lower && x > 0.) || (upper && x <= 0.)) \
202 *ccum = log1p(-exp(-xsq * xsq * 0.5) * \
203 exp(-del * 0.5) * temp); \
204 } \
205 else { \
206 *cum = exp(-xsq * xsq * 0.5) * exp(-del * 0.5) * temp; \
207 *ccum = 1.0 - *cum; \
208 }
209
210 #define swap_tail \
211 if (x > 0.) {/* swap ccum <--> cum */ \
212 temp = *cum; if(lower) *cum = *ccum; *ccum = temp; \
213 }
214
215 do_del(y);
216 swap_tail;
217 }
218
219 /* else |x| > sqrt(32) = 5.657 :
220 * the next two case differentiations were really for lower=T, log=F
221 * Particularly *not* for log_p !
222
223 * Cody had (-37.5193 < x && x < 8.2924) ; R originally had y < 50
224 *
225 * Note that we do want symmetry(0), lower/upper -> hence use y
226 */
227 else if((log_p && y < 1e170) /* avoid underflow below */
228 /* ^^^^^ MM FIXME: can speedup for log_p and much larger |x| !
229 * Then, make use of Abramowitz & Stegun, 26.2.13, something like
230
231 xsq = x*x;
232
233 if(xsq * DBL_EPSILON < 1.)
234 del = (1. - (1. - 5./(xsq+6.)) / (xsq+4.)) / (xsq+2.);
235 else
236 del = 0.;
237 *cum = -.5*xsq - M_LN_SQRT_2PI - log(x) + log1p(-del);
238 *ccum = log1p(-exp(*cum)); /.* ~ log(1) = 0 *./
239
240 swap_tail;
241
242 [Yes, but xsq might be infinite.]
243
244 */
245 || (lower && -37.5193 < x && x < 8.2924)
246 || (upper && -8.2924 < x && x < 37.5193)
247 ) {
248
249 /* Evaluate pnorm for x in (-37.5, -5.657) union (5.657, 37.5) */
250 xsq = 1.0 / (x * x); /* (1./x)*(1./x) might be better */
251 xnum = p[5] * xsq;
252 xden = xsq;
253 for (i = 0; i < 4; ++i) {
254 xnum = (xnum + p[i]) * xsq;
255 xden = (xden + q[i]) * xsq;
256 }
257 temp = xsq * (xnum + p[4]) / (xden + q[4]);
258 temp = (M_1_SQRT_2PI - temp) / y;
259
260 do_del(x);
261 swap_tail;
262 } else { /* large x such that probs are 0 or 1 */
263 if(x > 0) { *cum = R_D__1; *ccum = R_D__0; }
264 else { *cum = R_D__0; *ccum = R_D__1; }
265 }
266
267
268 #ifdef NO_DENORMS
269 /* do not return "denormalized" -- we do in R */
270 if(log_p) {
271 if(*cum > -min) *cum = -0.;
272 if(*ccum > -min)*ccum = -0.;
273 }
274 else {
275 if(*cum < min) *cum = 0.;
276 if(*ccum < min) *ccum = 0.;
277 }
278 #endif
279 return;
280 }
281