1 |
/* cdf/inverse_normal.c |
/* cdf/inverse_normal.c |
2 |
* |
* |
3 |
* Copyright (C) 2002 Przemyslaw Sliwa and Jason H. Stover. |
* Copyright (C) 2002 Przemyslaw Sliwa and Jason H. Stover. |
4 |
* |
* |
5 |
* This program is free software; you can redistribute it and/or modify |
* This program is free software; you can redistribute it and/or modify |
6 |
* it under the terms of the GNU General Public License as published by |
* it under the terms of the GNU General Public License as published by |
7 |
* the Free Software Foundation; either version 2 of the License, or (at |
* the Free Software Foundation; either version 2 of the License, or (at |
8 |
* your option) any later version. |
* your option) any later version. |
9 |
* |
* |
10 |
* This program is distributed in the hope that it will be useful, but |
* This program is distributed in the hope that it will be useful, but |
11 |
* WITHOUT ANY WARRANTY; without even the implied warranty of |
* WITHOUT ANY WARRANTY; without even the implied warranty of |
12 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
13 |
* General Public License for more details. |
* General Public License for more details. |
14 |
* |
* |
15 |
* You should have received a copy of the GNU General Public License |
* You should have received a copy of the GNU General Public License |
16 |
* along with this program; if not, write to the Free Software |
* along with this program; if not, write to the Free Software |
17 |
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA. |
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA. |
18 |
*/ |
*/ |
19 |
|
|
20 |
/* |
/* |
21 |
* Computes the inverse normal cumulative distribution function |
* Computes the inverse normal cumulative distribution function |
22 |
* according to the algorithm shown in |
* according to the algorithm shown in |
23 |
* |
* |
24 |
* Wichura, M.J. (1988). |
* Wichura, M.J. (1988). |
25 |
* Algorithm AS 241: The Percentage Points of the Normal Distribution. |
* Algorithm AS 241: The Percentage Points of the Normal Distribution. |
26 |
* Applied Statistics, 37, 477-484. |
* Applied Statistics, 37, 477-484. |
27 |
*/ |
*/ |
28 |
|
|
29 |
#include <config.h> |
#include <config.h> |
30 |
#include <gsl/gsl_errno.h> |
#include <gsl/gsl_errno.h> |
31 |
#include <gsl/gsl_math.h> |
#include <gsl/gsl_math.h> |
32 |
#include <gsl/gsl_cdf.h> |
#include <gsl/gsl_cdf.h> |
33 |
#include <gsl/gsl_cdf_result.h> |
|
34 |
#include "error.h" |
#include "rat_eval.h" |
35 |
|
|
36 |
int gsl_cdf_inv_normal_e(double p, gsl_cdf_result * result) |
double |
37 |
{ |
gsl_cdf_ugaussian_P_inv (double p) |
38 |
double q, r, ret_val; |
{ |
39 |
|
double r, x; |
40 |
q = p - .5; |
|
41 |
if (fabs (q) <= .425) |
double q = p - 0.5; |
42 |
{ |
|
43 |
r = .180625 - q * q; |
if (fabs (q) <= 0.425) |
44 |
ret_val = q * (((((((r * 2509.0809287301226727 + |
{ |
45 |
33430.575583588128105) * r + |
const double a[8] = { 3.387132872796366608, 133.14166789178437745, |
46 |
67265.770927008700853) * r + |
1971.5909503065514427, 13731.693765509461125, |
47 |
45921.953931549871457) * r + |
45921.953931549871457, 67265.770927008700853, |
48 |
13731.693765509461125) * r + |
33430.575583588128105, 2509.0809287301226727 |
49 |
1971.5909503065514427) * r + |
}; |
50 |
133.14166789178437745) * r + |
|
51 |
3.387132872796366608) / |
const double b[8] = { 1.0, 42.313330701600911252, |
52 |
(((((((r |
687.1870074920579083, 5394.1960214247511077, |
53 |
* 5226.495278852854561 + 28729.085735721942674) * r + |
21213.794301586595867, 39307.89580009271061, |
54 |
39307.89580009271061) * r + 21213.794301586595867) * r + |
28729.085735721942674, 5226.495278852854561 |
55 |
5394.1960214247511077) * r + 687.1870074920579083) * r + |
}; |
56 |
42.313330701600911252) * r + 1.); |
|
57 |
result->val = ret_val; |
r = 0.180625 - q * q; |
58 |
return GSL_SUCCESS; |
|
59 |
} |
x = q * rat_eval (a, 8, b, 8, r); |
60 |
|
|
61 |
else |
return x; |
62 |
{ |
} |
63 |
if (q < 0.) |
|
64 |
{ |
if (q < 0.0) |
65 |
r = p; |
{ |
66 |
} |
r = p; |
67 |
else |
} |
68 |
{ |
else |
69 |
r = 1.0 - p; |
{ |
70 |
} |
r = 1.0 - p; |
71 |
if (r <= 0.0) |
} |
72 |
{ |
|
73 |
DOMAIN_ERROR(result); |
r = sqrt (-log (r)); |
74 |
} |
|
75 |
r = sqrt (-log(r)); |
if (r <= 5.0) |
76 |
if (r <= 5.0) |
{ |
77 |
{ |
const double a[] = { 1.42343711074968357734, 4.6303378461565452959, |
78 |
r += -1.6; |
5.7694972214606914055, 3.64784832476320460504, |
79 |
ret_val = (((((((r * 7.7454501427834140764e-4 + |
1.27045825245236838258, 0.24178072517745061177, |
80 |
.0227238449892691845833) * r + |
0.0227238449892691845833, 7.7454501427834140764e-4 |
81 |
.24178072517745061177) * r + |
}; |
82 |
1.27045825245236838258) * r + |
|
83 |
3.64784832476320460504) * r + |
const double b[] = { 1.0, 2.05319162663775882187, |
84 |
5.7694972214606914055) * r + |
1.6763848301838038494, 0.68976733498510000455, |
85 |
4.6303378461565452959) * r + |
0.14810397642748007459, 0.0151986665636164571966, |
86 |
1.42343711074968357734) / |
5.475938084995344946e-4, 1.05075007164441684324e-9 |
87 |
(((((((r |
}; |
88 |
* 1.05075007164441684324e-9 + |
|
89 |
5.475938084995344946e-4) * r + |
x = rat_eval (a, 8, b, 8, (r - 1.6)); |
90 |
.0151986665636164571966) * r + |
} |
91 |
.14810397642748007459) * r + .68976733498510000455) * r + |
else |
92 |
1.6763848301838038494) * r + 2.05319162663775882187) * r + |
{ |
93 |
1.0); |
const double a[] = { 6.6579046435011037772, 5.4637849111641143699, |
94 |
} |
1.7848265399172913358, 0.29656057182850489123, |
95 |
else |
0.026532189526576123093, 0.0012426609473880784386, |
96 |
{ |
2.71155556874348757815e-5, 2.01033439929228813265e-7 |
97 |
r += -5.0; |
}; |
98 |
ret_val = (((((((r * 2.01033439929228813265e-7 + |
|
99 |
2.71155556874348757815e-5) * r + |
const double b[] = { 1.0, 0.59983220655588793769, |
100 |
.0012426609473880784386) * r + |
0.13692988092273580531, 0.0148753612908506148525, |
101 |
.026532189526576123093) * r + |
7.868691311456132591e-4, 1.8463183175100546818e-5, |
102 |
.29656057182850489123) * r + |
1.4215117583164458887e-7, 2.04426310338993978564e-15 |
103 |
1.7848265399172913358) * r + |
}; |
104 |
5.4637849111641143699) * r + |
|
105 |
6.6579046435011037772) / |
x = rat_eval (a, 8, b, 8, (r - 5.0)); |
106 |
(((((((r |
} |
107 |
* 2.04426310338993978564e-15 + |
|
108 |
1.4215117583164458887e-7) * r + |
if (q < 0.0) |
109 |
1.8463183175100546818e-5) * r + |
{ |
110 |
7.868691311456132591e-4) * r + |
return -x; |
111 |
.0148753612908506148525) * r + |
} |
112 |
.13692988092273580531) * r + .59983220655588793769) * r + |
else |
113 |
1.0); |
{ |
114 |
} |
return x; |
115 |
if (q < 0.0) |
} |
116 |
{ |
|
117 |
ret_val = -ret_val; |
} |
118 |
} |
|
119 |
result->val = ret_val; |
double |
120 |
return GSL_SUCCESS; |
gsl_cdf_ugaussian_Q_inv (double p) |
121 |
} |
{ |
122 |
result->val = ret_val; |
return gsl_cdf_ugaussian_P_inv (1.0 - p); |
123 |
return GSL_SUCCESS; |
} |
|
|
|
|
} |
|
|
/* |
|
|
* No error code. |
|
|
*/ |
|
|
#include "eval.h" |
|
|
|
|
|
double gsl_cdf_inv_normal ( double p ) |
|
|
{ |
|
|
EVAL_RESULT ( gsl_cdf_inv_normal_e ( p, &result)); |
|
|
} |
|
|
|
|
|
/* int main(void) */ |
|
|
|
|
|
/* { */ |
|
|
|
|
|
/* double quantile, p; */ |
|
|
/* int ifail; */ |
|
|
|
|
|
/* for(p=1.0E-15; p<1.0; p+=1.0E-5) */ |
|
|
/* { */ |
|
|
/* quantile = gsl_normal_quantile(p, &ifail); */ |
|
|
/* printf("VALUE: %.19f\tIFAIL %d\n", quantile, ifail); */ |
|
|
/* } */ |
|
|
|
|
|
/* return 0; */ |
|
|
|
|
|
/* } */ |
|