36 GaussianDistribution::GaussianDistribution(
const double mean,
const double standardDeviation) {
38 p_stdev = standardDeviation;
49 double GaussianDistribution::Probability(
const double value) {
51 return std::exp(-0.5 * std::pow((value - p_mean) / p_stdev, 2)) / (std::sqrt(2 *
PI) * p_stdev);
62 double GaussianDistribution::CumulativeDistribution(
const double value) {
64 if(value == DBL_MIN) {
67 else if(value == DBL_MAX) {
73 double x = (value - p_mean) / p_stdev;
80 for(
int n = 0; pre != sum; n++) {
83 sum += std::pow(x, 2 * n + 1) / (fact * (2 * n + 1) * std::pow(-2.0, n));
88 return 50.0 + 100.0 / std::sqrt(2.0 *
PI) * sum;
104 double GaussianDistribution::InverseCumulativeDistribution(
const double percent) {
106 static double lowCutoff = 2.425;
107 static double highCutoff = 97.575;
109 if((percent < 0.0) || (percent > 100.0)) {
110 string m =
"Argument percent outside of the range 0 to 100 in";
111 m +=
" [GaussianDistribution::InverseCumulativeDistribution]";
123 else if(percent == 100.0) {
127 if(percent < lowCutoff) {
128 double q = std::sqrt(-2.0 * std::log(percent / 100.0));
131 else if(percent < highCutoff) {
132 double q = percent / 100.0 - 0.5,
137 double q = std::sqrt(-2.0 * std::log(1.0 - percent / 100.0));
138 x = -1.0 * C(q) / D(q);
143 double e = CumulativeDistribution(p_stdev * x + p_mean) - percent;
144 double u = e * std::sqrt(2.0 *
PI) * std::exp(-0.5 * x * x);
145 x = x - u / (1 + 0.5 * x * u);
147 return p_stdev * x + p_mean;
154 double GaussianDistribution::A(
const double x) {
155 static const double a[6] = { -39.69683028665376,
163 return((((a[0] * x + a[1]) * x + a[2]) * x + a[3]) * x + a[4]) * x + a[5];
166 double GaussianDistribution::B(
const double x) {
167 static const double b[6] = { -54.47609879822406,
175 return((((b[0] * x + b[1]) * x + b[2]) * x + b[3]) * x + b[4]) * x + b[5];
178 double GaussianDistribution::C(
const double x) {
179 static const double c[6] = { -0.007784894002430293,
187 return((((c[0] * x + c[1]) * x + c[2]) * x + c[3]) * x + c[4]) * x + c[5];
190 double GaussianDistribution::D(
const double x) {
191 static const double d[5] = { 0.007784695709041462,
198 return(((d[0] * x + d[1]) * x + d[2]) * x + d[3]) * x + d[4];