64 lines
1.1 KiB
C
64 lines
1.1 KiB
C
#include <math.h>
|
|
#include <stdlib.h>
|
|
#include <time.h>
|
|
|
|
|
|
#define M_PI 3.14159265358979323846264338327950288
|
|
|
|
|
|
double normalcdf(double z)
|
|
{
|
|
/* https://www.johndcook.com/blog/cpp_phi/ */
|
|
double a1, a2, a3, a4, a5, p, t, y;
|
|
int sign;
|
|
|
|
a1 = 0.254829592;
|
|
a2 = -0.284496736;
|
|
a3 = 1.421413741;
|
|
a4 = -1.453152027;
|
|
a5 = 1.061405429;
|
|
p = 0.3275911;
|
|
|
|
sign = 1;
|
|
if (z < 0)
|
|
sign = -1;
|
|
z = fabs(z) / pow(2, 0.5);
|
|
t = 1 / (1 + p * z);
|
|
y = 1 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * exp(-z * z);
|
|
return 0.5 * (1 + sign * y);
|
|
}
|
|
|
|
|
|
double normalpdf(double z)
|
|
{
|
|
return exp(-0.5 * pow(z, 2)) / pow(M_PI * 2, 0.5);
|
|
}
|
|
|
|
|
|
double gaussrand(void)
|
|
{
|
|
/* Marsaglia and Bray, ``A Convenient Method for Generating Normal Variables'' */
|
|
static double V1, V2, S;
|
|
static int phase;
|
|
double X;
|
|
|
|
if (phase == 0) {
|
|
do {
|
|
double U1 = (double)rand() / RAND_MAX;
|
|
double U2 = (double)rand() / RAND_MAX;
|
|
|
|
V1 = 2 * U1 - 1;
|
|
V2 = 2 * U2 - 1;
|
|
S = V1 * V1 + V2 * V2;
|
|
} while (S >= 1 || S == 0);
|
|
|
|
X = V1 * sqrt(-2 * log(S) / S);
|
|
} else {
|
|
X = V2 * sqrt(-2 * log(S) / S);
|
|
}
|
|
|
|
phase = 1 - phase;
|
|
|
|
return X;
|
|
}
|