NormalDistribution.cc

Go to the documentation of this file.

#include <cmath>
#include <limits>
#include <utl/NormalDistribution.h>
#include <utl/ErrorLogger.h>

using namespace std;


namespace utl {

  double
  InverseNormalCDF(const double p)
  {
    const double a[6] = {
       -39.69683028665376,
       220.9460984245205,
      -275.9285104469687,
       138.3577518672690,
       -30.66479806614716,
         2.506628277459239
    };

    const double b[5] = {
       -54.47609879822406,
       161.5858368580409,
      -155.6989798598866,
        66.80131188771972,
       -13.28068155288572
    };

    const double c[6] = {
      -0.007784894002430293,
      -0.3223964580411365,
      -2.400758277161838,
      -2.549732539343734,
       4.374664141464968,
       2.938163982698783
    };

    const double d[4] = {
      0.007784695709041462,
      0.3224671290700398,
      2.445134137142996,
      3.754408661907416
    };

    if (p < 0 || p > 1) {
      ERROR("Bad p in routine InverseNormalCDF");
      return numeric_limits<double>::quiet_NaN();
    }

    if (p == 0)
      return -numeric_limits<double>::infinity();

    if (p == 1)
      return numeric_limits<double>::infinity();

    // inverse normal cdf is an odd function

    const double pp = min(p, 1 - p);

    double x;

    // minmax rational approximations are better than 8 decimal places    
    if (pp <= 0.02425) {
      const double q = sqrt(-2 * log(pp));
      x =
        (((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5]) /
        ((((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1);
    } else {
      const double q = pp - 0.5;
      const double r = q*q;
      x =
        (((((a[0]*r + a[1])*r + a[2])*r + a[3])*r + a[4])*r + a[5])*q /
        (((((b[0]*r + b[1])*r + b[2])*r + b[3])*r + b[4])*r + 1);
    }

    // one step of Halley's iteration brings accuracy to the full machine precision.
    const double df = NormalCDF(x) - pp;
    const double dfof = df * sqrt(2*kPi) * exp(0.5*Sqr(x));
    x -= dfof / (1 + 0.5*x*dfof);

    return p > 0.5 ? -x : x;
  }


  double
  TruncatedNormalPDF(const double x, const double mu, const double sigma,
                     const double xmin, const double xmax)
  {
    if (xmin >= xmax || x < xmin || x > xmax)
      return 0;
    return NormalPDF(x, mu, sigma) / (NormalCDF(xmax, mu, sigma) - NormalCDF(xmin, mu, sigma));
  }


  // truncated gaussian with xmin = 0 and xmax = infinity
  double
  TruncatedNormalPDF(const double x, const double mu, const double sigma, const double xmin)
  {
    if (x <= xmin)
      return 0;
    return NormalPDF(x, mu, sigma) / (1 - NormalCDF(xmin, mu, sigma));
  }

}