1 /*
2 * Mathlib : A C Library of Special Functions
3 * Copyright (C) 2000--2020 The R Core Team
4 * Copyright (C) 1998 Ross Ihaka
5 * based on AS 111 (C) 1977 Royal Statistical Society
6 * and on AS 241 (C) 1988 Royal Statistical Society
7 *
8 * This program is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with this program; if not, a copy is available at
20 * https://www.R-project.org/Licenses/
21 *
22 * SYNOPSIS
23 *
24 * double qnorm5(double p, double mu, double sigma,
25 * int lower_tail, int log_p)
26 * {qnorm (..) is synonymous and preferred inside R}
27 *
28 * DESCRIPTION
29 *
30 * Compute the quantile function for the normal distribution.
31 *
32 * For small to moderate probabilities, algorithm referenced
33 * below is used to obtain an initial approximation which is
34 * polished with a final Newton step.
35 *
36 * For very large arguments, an algorithm of Wichura is used.
37 *
38 * REFERENCE
39 *
40 * Beasley, J. D. and S. G. Springer (1977).
41 * Algorithm AS 111: The percentage points of the normal distribution,
42 * Applied Statistics, 26, 118-121.
43 *
44 * Wichura, M.J. (1988).
45 * Algorithm AS 241: The Percentage Points of the Normal Distribution.
46 * Applied Statistics, 37, 477-484.
47 */
48
49 #include "nmath.h"
50 #include "dpq.h"
51
qnorm5(double p,double mu,double sigma,int lower_tail,int log_p)52 double qnorm5(double p, double mu, double sigma, int lower_tail, int log_p)
53 {
54 double p_, q, r, val;
55
56 #ifdef IEEE_754
57 if (ISNAN(p) || ISNAN(mu) || ISNAN(sigma))
58 return p + mu + sigma;
59 #endif
60 R_Q_P01_boundaries(p, ML_NEGINF, ML_POSINF);
61
62 if(sigma < 0) ML_WARN_return_NAN;
63 if(sigma == 0) return mu;
64
65 p_ = R_DT_qIv(p);/* real lower_tail prob. p */
66 q = p_ - 0.5;
67
68 #ifdef DEBUG_qnorm
69 REprintf("qnorm(p=%10.7g, m=%g, s=%g, l.t.= %d, log= %d): q = %g\n",
70 p,mu,sigma, lower_tail, log_p, q);
71 #endif
72
73
74 /*-- use AS 241 --- */
75 /* double ppnd16_(double *p, long *ifault)*/
76 /* ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3
77
78 Produces the normal deviate Z corresponding to a given lower
79 tail area of P; Z is accurate to about 1 part in 10**16.
80
81 (original fortran code used PARAMETER(..) for the coefficients
82 and provided hash codes for checking them...)
83 */
84 if (fabs(q) <= .425) {/* |p~ - 0.5| <= .425 <==> 0.075 <= p~ <= 0.925 */
85 r = .180625 - q * q; // = .425^2 - q^2 >= 0
86 val =
87 q * (((((((r * 2509.0809287301226727 +
88 33430.575583588128105) * r + 67265.770927008700853) * r +
89 45921.953931549871457) * r + 13731.693765509461125) * r +
90 1971.5909503065514427) * r + 133.14166789178437745) * r +
91 3.387132872796366608)
92 / (((((((r * 5226.495278852854561 +
93 28729.085735721942674) * r + 39307.89580009271061) * r +
94 21213.794301586595867) * r + 5394.1960214247511077) * r +
95 687.1870074920579083) * r + 42.313330701600911252) * r + 1.);
96 }
97 else { /* closer than 0.075 from {0,1} boundary :
98 * r := log(p~); p~ = min(p, 1-p) < 0.075 : */
99 if(log_p && ((lower_tail && q <= 0) || (!lower_tail && q > 0))) {
100 r = p;
101 } else {
102 r = log( (q > 0) ? R_DT_CIv(p) /* 1-p */ : p_ /* = R_DT_Iv(p) ^= p */);
103 }
104 // r = sqrt( - log(min(p,1-p)) ) <==> min(p, 1-p) = exp( - r^2 ) :
105 r = sqrt(-r);
106 #ifdef DEBUG_qnorm
107 REprintf("\t close to 0 or 1: r = %7g\n", r);
108 #endif
109 if (r <= 5.) { /* <==> min(p,1-p) >= exp(-25) ~= 1.3888e-11 */
110 r += -1.6;
111 val = (((((((r * 7.7454501427834140764e-4 +
112 .0227238449892691845833) * r + .24178072517745061177) *
113 r + 1.27045825245236838258) * r +
114 3.64784832476320460504) * r + 5.7694972214606914055) *
115 r + 4.6303378461565452959) * r +
116 1.42343711074968357734)
117 / (((((((r *
118 1.05075007164441684324e-9 + 5.475938084995344946e-4) *
119 r + .0151986665636164571966) * r +
120 .14810397642748007459) * r + .68976733498510000455) *
121 r + 1.6763848301838038494) * r +
122 2.05319162663775882187) * r + 1.);
123 }
124 else if(r >= 816) { // p is *extremly* close to 0 or 1 - only possibly when log_p =TRUE
125 // Using the asymptotical formula -- is *not* optimal but uniformly better than branch below
126 val = r * M_SQRT2;
127 }
128 else { // p is very close to 0 or 1: r > 5 <==> min(p,1-p) < exp(-25) = 1.3888..e-11
129 r += -5.;
130 val = (((((((r * 2.01033439929228813265e-7 +
131 2.71155556874348757815e-5) * r +
132 .0012426609473880784386) * r + .026532189526576123093) *
133 r + .29656057182850489123) * r +
134 1.7848265399172913358) * r + 5.4637849111641143699) *
135 r + 6.6579046435011037772)
136 / (((((((r *
137 2.04426310338993978564e-15 + 1.4215117583164458887e-7)*
138 r + 1.8463183175100546818e-5) * r +
139 7.868691311456132591e-4) * r + .0148753612908506148525)
140 * r + .13692988092273580531) * r +
141 .59983220655588793769) * r + 1.);
142 }
143
144 if(q < 0.0)
145 val = -val;
146 /* return (q >= 0.)? r : -r ;*/
147 }
148 return mu + sigma * val;
149 }
150