Loading...
Searching...
No Matches
normp.c
Go to the documentation of this file.
#include <math.h>
/*-
* SUBROUTINE Cdhc_normp(Z, P, Q, PDF)
*
* Normal distribution probabilities accurate to 1.e-15.
* Z = no. of standard deviations from the mean.
* P, Q = probabilities to the left & right of Z. P + Q = 1.
* PDF = the probability density.
*
* Based upon algorithm 5666 for the error function, from:
* Hart, J.F. et al, 'Computer Approximations', Wiley 1968
*
* Programmer: Alan Miller
*
* Latest revision - 30 March 1986
*
*/
/* Conversion to C by James Darrell McCauley, 24 September 1994 */
double Cdhc_normp(double z)
{
static double p[7] = {220.2068679123761, 221.2135961699311,
112.079291497870, 33.91286607838300,
6.37396220353165, 0.7003830644436881,
0.352624965998910e-1};
static double q[8] = {440.4137358247522, 793.8265125199484,
637.3336333788311, 296.5642487796737,
86.78073220294608, 16.06417757920695,
1.755667163182642, 0.8838834764831844e-1};
static double cutoff = 7.071, root2pi = 2.506628274631001;
double zabs, expntl;
double pee, queue, pdf;
zabs = fabs(z);
if (zabs > 37.0) {
pdf = 0.0;
if (z > 0.0) {
pee = 1.0;
queue = 0.0;
}
else {
pee = 0.0;
queue = 1.0;
}
return pee;
}
expntl = exp(-0.5 * zabs * zabs);
pdf = expntl / root2pi;
if (zabs < cutoff)
pee = expntl *
((((((p[6] * zabs + p[5]) * zabs + p[4]) * zabs + p[3]) * zabs +
p[2]) *
zabs +
p[1]) *
zabs +
p[0]) /
(((((((q[7] * zabs + q[6]) * zabs + q[5]) * zabs + q[4]) * zabs +
q[3]) *
zabs +
q[2]) *
zabs +
q[1]) *
zabs +
q[0]);
else
pee =
pdf /
(zabs +
1.0 / (zabs + 2.0 / (zabs + 3.0 / (zabs + 4.0 / (zabs + 0.65)))));
if (z < 0.0)
queue = 1.0 - pee;
else {
queue = pee;
pee = 1.0 - queue;
}
return pee;
}
Cdhc_normp
double Cdhc_normp(double z)
Definition normp.c:22