46 |
#include <gsl/gsl_math.h> |
#include <gsl/gsl_math.h> |
47 |
#include <gsl/gsl_cdf.h> |
#include <gsl/gsl_cdf.h> |
48 |
|
|
49 |
#ifndef M_1_SQRT2PI |
#ifndef M_1_SQRT2PI |
50 |
#define M_1_SQRT2PI (M_2_SQRTPI * M_SQRT1_2 / 2.0) |
#define M_1_SQRT2PI (M_2_SQRTPI * M_SQRT1_2 / 2.0) |
51 |
#endif |
#endif |
52 |
|
|
53 |
|
#define SQRT32 (4.0 * M_SQRT2) |
54 |
|
|
55 |
/* |
/* |
56 |
* IEEE double precision dependent constants. |
* IEEE double precision dependent constants. |
57 |
* |
* |
62 |
*/ |
*/ |
63 |
|
|
64 |
#define GAUSS_EPSILON GSL_DBL_EPSILON / 2 |
#define GAUSS_EPSILON GSL_DBL_EPSILON / 2 |
65 |
#define GAUSS_XUPPER 8.572 |
#define GAUSS_XUPPER 8.572 |
66 |
#define GAUSS_XLOWER -37.519 |
#define GAUSS_XLOWER -37.519 |
67 |
|
|
68 |
#define GAUSS_SCALE 16.0 |
#define GAUSS_SCALE 16.0 |
69 |
|
|
70 |
static double get_del ( double x, double rational ) |
static double |
71 |
|
get_del (double x, double rational) |
72 |
{ |
{ |
73 |
double xsq = 0.0; |
double xsq = 0.0; |
74 |
double del = 0.0; |
double del = 0.0; |
75 |
double result = 0.0; |
double result = 0.0; |
76 |
|
|
77 |
xsq = floor( x * GAUSS_SCALE ) / GAUSS_SCALE; |
xsq = floor (x * GAUSS_SCALE) / GAUSS_SCALE; |
78 |
del = (x - xsq) * (x + xsq); |
del = (x - xsq) * (x + xsq); |
79 |
del *= 0.5; |
del *= 0.5; |
80 |
|
|
81 |
result = exp(-0.5 * xsq * xsq ) * exp(-1.0 * del) * rational; |
result = exp (-0.5 * xsq * xsq) * exp (-1.0 * del) * rational; |
82 |
|
|
83 |
return result; |
return result; |
84 |
} |
} |
85 |
|
|
86 |
/* |
/* |
87 |
* Normal cdf for fabs(x) < 0.66291 |
* Normal cdf for fabs(x) < 0.66291 |
88 |
*/ |
*/ |
89 |
static double gauss_smallx(const double x) |
static double |
90 |
|
gauss_smallx (const double x) |
91 |
{ |
{ |
92 |
|
unsigned int i; |
93 |
double result = 0.0; |
double result = 0.0; |
|
double temp = 0.0; |
|
94 |
double xsq; |
double xsq; |
95 |
double xnum; |
double xnum; |
96 |
double xden; |
double xden; |
108 |
10260.932208618978205, |
10260.932208618978205, |
109 |
45507.789335026729956 |
45507.789335026729956 |
110 |
}; |
}; |
|
int n; |
|
111 |
|
|
112 |
xsq = x * x; |
xsq = x * x; |
113 |
xnum = a[4] * xsq; |
xnum = a[4] * xsq; |
114 |
xden = xsq; |
xden = xsq; |
115 |
for (n = 0; n < 3; ++n) { |
|
116 |
xnum = (xnum + a[n]) * xsq; |
for (i = 0; i < 3; i++) |
117 |
xden = (xden + b[n]) * xsq; |
{ |
118 |
} |
xnum = (xnum + a[i]) * xsq; |
119 |
|
xden = (xden + b[i]) * xsq; |
120 |
|
} |
121 |
|
|
122 |
result = x * (xnum + a[3]) / (xden + b[3]); |
result = x * (xnum + a[3]) / (xden + b[3]); |
|
|
|
|
return result; |
|
123 |
|
|
124 |
|
return result; |
125 |
} |
} |
126 |
|
|
127 |
/* |
/* |
128 |
* Normal cdf for 0.66291 < fabs(x) < sqrt(32). |
* Normal cdf for 0.66291 < fabs(x) < sqrt(32). |
129 |
*/ |
*/ |
130 |
static double gauss_mediumx ( const double x) |
static double |
131 |
|
gauss_mediumx (const double x) |
132 |
{ |
{ |
133 |
int n; |
unsigned int i; |
134 |
double temp = 0.0; |
double temp = 0.0; |
135 |
double result = 0.0; |
double result = 0.0; |
136 |
double xnum; |
double xnum; |
158 |
38912.003286093271411, |
38912.003286093271411, |
159 |
19685.429676859990727 |
19685.429676859990727 |
160 |
}; |
}; |
161 |
absx = fabs(x); |
|
162 |
|
absx = fabs (x); |
163 |
|
|
164 |
xnum = c[8] * absx; |
xnum = c[8] * absx; |
165 |
xden = absx; |
xden = absx; |
166 |
for (n = 0; n < 7; ++n) { |
|
167 |
xnum = (xnum + c[n]) * absx; |
for (i = 0; i < 7; i++) |
168 |
xden = (xden + d[n]) * absx; |
{ |
169 |
} |
xnum = (xnum + c[i]) * absx; |
170 |
|
xden = (xden + d[i]) * absx; |
171 |
|
} |
172 |
|
|
173 |
temp = (xnum + c[7]) / (xden + d[7]); |
temp = (xnum + c[7]) / (xden + d[7]); |
174 |
|
|
175 |
result = get_del(x, temp); |
result = get_del (x, temp); |
176 |
|
|
177 |
return result; |
return result; |
178 |
} |
} |
181 |
* Normal cdf for |
* Normal cdf for |
182 |
* {sqrt(32) < x < GAUSS_XUPPER} union { GAUSS_XLOWER < x < -sqrt(32) }. |
* {sqrt(32) < x < GAUSS_XUPPER} union { GAUSS_XLOWER < x < -sqrt(32) }. |
183 |
*/ |
*/ |
184 |
static double gauss_grandex ( const double x ) |
static double |
185 |
|
gauss_grandex (const double x) |
186 |
{ |
{ |
187 |
int n; |
int i; |
188 |
double result; |
double result; |
189 |
double xsq; |
double xsq; |
190 |
double temp; |
double temp; |
191 |
double xnum; |
double xnum; |
192 |
double xden; |
double xden; |
193 |
double absx; |
double absx; |
194 |
|
|
195 |
const double p[6] = { |
const double p[6] = { |
196 |
0.21589853405795699, |
0.21589853405795699, |
197 |
0.1274011611602473639, |
0.1274011611602473639, |
208 |
7.29751555083966205e-5 |
7.29751555083966205e-5 |
209 |
}; |
}; |
210 |
|
|
211 |
absx = fabs(x); |
absx = fabs (x); |
212 |
xsq = 1.0 / (x * x); |
xsq = 1.0 / (x * x); |
213 |
xnum = p[5] * xsq; |
xnum = p[5] * xsq; |
214 |
xden = xsq; |
xden = xsq; |
215 |
for (n = 0; n < 4; ++n) { |
|
216 |
xnum = (xnum + p[n]) * xsq; |
for (i = 0; i < 4; i++) |
217 |
xden = (xden + q[n]) * xsq; |
{ |
218 |
} |
xnum = (xnum + p[i]) * xsq; |
219 |
|
xden = (xden + q[i]) * xsq; |
220 |
|
} |
221 |
|
|
222 |
temp = xsq * (xnum + p[4]) / (xden + q[4]); |
temp = xsq * (xnum + p[4]) / (xden + q[4]); |
223 |
temp = (M_1_SQRT2PI - temp) / absx; |
temp = (M_1_SQRT2PI - temp) / absx; |
224 |
|
|
225 |
result = get_del(x, temp); |
result = get_del (x, temp); |
226 |
|
|
227 |
return result; |
return result; |
228 |
} |
} |
229 |
|
|
230 |
double gsl_cdf_gaussian_Q (const double x ) |
double |
231 |
|
gsl_cdf_gaussian_P (const double x) |
232 |
{ |
{ |
233 |
double result; |
double result; |
234 |
double absx; |
double absx = fabs (x); |
|
double sqrt32 = sqrt(32); |
|
|
|
|
|
absx = fabs(x); |
|
235 |
|
|
236 |
if ( absx < 0.66291 ) |
if (absx < GAUSS_EPSILON) |
237 |
{ |
{ |
238 |
result = gauss_smallx( x); |
result = 0.5; |
239 |
if ( x < 0.0 ) |
return result; |
|
{ |
|
|
result = fabs(result) + 0.5; |
|
|
} |
|
|
else |
|
|
result = 0.5 - result; |
|
|
|
|
240 |
} |
} |
241 |
else if ( absx < sqrt(32)) |
else if (absx < 0.66291) |
242 |
{ |
{ |
243 |
result = gauss_mediumx(x); |
result = 0.5 + gauss_smallx (x); |
244 |
if( x < 0.0 ) |
return result; |
|
{ |
|
|
result = 1.0 - result; |
|
|
} |
|
245 |
} |
} |
246 |
else if ( x >= -GAUSS_XUPPER || x <= -GAUSS_XLOWER ) |
else if (absx < SQRT32) |
247 |
{ |
{ |
248 |
result = gauss_grandex (x); |
result = gauss_mediumx (x); |
249 |
if ( x < 0.0 ) |
|
250 |
|
if (x > 0.0) |
251 |
{ |
{ |
252 |
result = 1.0 - result; |
result = 1.0 - result; |
253 |
} |
} |
254 |
|
|
255 |
|
return result; |
256 |
} |
} |
257 |
else if ( x > -GAUSS_XLOWER ) |
else if (x > GAUSS_XUPPER) |
258 |
{ |
{ |
259 |
result = 0.0; |
result = 1.0; |
260 |
|
return result; |
261 |
} |
} |
262 |
else if ( x < -GAUSS_XUPPER ) |
else if (x < GAUSS_XLOWER) |
263 |
{ |
{ |
264 |
result = 1.0; |
result = 0.0; |
265 |
|
return result; |
266 |
} |
} |
267 |
else if ( x > 0.0 && x < GAUSS_EPSILON ) |
else |
268 |
{ |
{ |
269 |
result = 0.5; |
result = gauss_grandex (x); |
270 |
|
|
271 |
|
if (x > 0.0) |
272 |
|
{ |
273 |
|
result = 1.0 - result; |
274 |
|
} |
275 |
} |
} |
276 |
|
|
277 |
return result; |
return result; |
278 |
} |
} |
279 |
double gsl_cdf_gaussian_P(const double x) |
|
280 |
|
double |
281 |
|
gsl_cdf_gaussian_Q (const double x) |
282 |
{ |
{ |
283 |
double result; |
double result; |
284 |
double absx; |
double absx = fabs (x); |
|
double sqrt32 = sqrt(32); |
|
285 |
|
|
286 |
absx = fabs(x); |
if (absx < GAUSS_EPSILON) |
|
|
|
|
if ( absx < 0.66291 ) |
|
287 |
{ |
{ |
288 |
result = 0.5 + gauss_smallx( x); |
result = 0.5; |
289 |
|
return result; |
290 |
} |
} |
291 |
else if ( absx < sqrt(32)) |
else if (absx < 0.66291) |
292 |
{ |
{ |
293 |
result = gauss_mediumx(x); |
result = gauss_smallx (x); |
294 |
if ( x > 0.0 ) |
|
295 |
|
if (x < 0.0) |
296 |
{ |
{ |
297 |
result = 1.0 - result; |
result = fabs (result) + 0.5; |
298 |
|
} |
299 |
|
else |
300 |
|
{ |
301 |
|
result = 0.5 - result; |
302 |
} |
} |
303 |
|
|
304 |
|
return result; |
305 |
} |
} |
306 |
else if ( x < GAUSS_XUPPER || x < GAUSS_XLOWER ) |
else if (absx < SQRT32) |
307 |
{ |
{ |
308 |
result = gauss_grandex (x); |
result = gauss_mediumx (x); |
309 |
if( x > 0.0 ) |
|
310 |
|
if (x < 0.0) |
311 |
{ |
{ |
312 |
result = 1.0 - result; |
result = 1.0 - result; |
313 |
} |
} |
314 |
|
|
315 |
|
return result; |
316 |
} |
} |
317 |
else if ( x > GAUSS_XUPPER ) |
else if (x > (-GAUSS_XLOWER)) |
318 |
{ |
{ |
319 |
result = 1.0; |
result = 0.0; |
320 |
|
return result; |
321 |
} |
} |
322 |
else if ( x < GAUSS_XLOWER ) |
else if (x < (-GAUSS_XUPPER)) |
323 |
{ |
{ |
324 |
result = 0.0; |
result = 1.0; |
325 |
|
return result; |
326 |
} |
} |
327 |
else if ( x > 0.0 && x < GAUSS_EPSILON ) |
else |
328 |
{ |
{ |
329 |
result = 0.5; |
result = gauss_grandex (x); |
330 |
|
|
331 |
|
if (x < 0.0) |
332 |
|
{ |
333 |
|
result = 1.0 - result; |
334 |
|
} |
335 |
|
|
336 |
} |
} |
337 |
|
|
338 |
return result; |
return result; |