inv_Phi.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_SCAL_FUN_INV_PHI_HPP
2 #define STAN_MATH_PRIM_SCAL_FUN_INV_PHI_HPP
3 
8 
9 namespace stan {
10 namespace math {
11 
12 /**
13  * The inverse of the unit normal cumulative distribution function.
14  *
15  * The return value for a specified input probability, $p$, is
16  * the unit normal variate, $x$, such that
17  *
18  * \f$\Phi(x) = \int_{-\infty}^x \mbox{\sf Norm}(x|0, 1) \ dx = p\f$
19  *
20  * Algorithm first derived in 2003 by Peter Jon Aklam at
21  * http://home.online.no/~pjacklam/notes/invnorm/
22  *
23  * @param p Argument between 0 and 1.
24  * @return Real number
25  */
26 inline double inv_Phi(double p) {
27  check_bounded("inv_Phi", "Probability variable", p, 0, 1);
28 
29  if (p < 8e-311)
30  return NEGATIVE_INFTY;
31  if (p == 1)
32  return INFTY;
33 
34  static const double a[6]
35  = {-3.969683028665376e+01, 2.209460984245205e+02, -2.759285104469687e+02,
36  1.383577518672690e+02, -3.066479806614716e+01, 2.506628277459239e+00};
37  static const double b[5]
38  = {-5.447609879822406e+01, 1.615858368580409e+02, -1.556989798598866e+02,
39  6.680131188771972e+01, -1.328068155288572e+01};
40  static const double c[6]
41  = {-7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e+00,
42  -2.549732539343734e+00, 4.374664141464968e+00, 2.938163982698783e+00};
43  static const double d[4] = {7.784695709041462e-03, 3.224671290700398e-01,
44  2.445134137142996e+00, 3.754408661907416e+00};
45 
46  static const double p_low = 0.02425;
47  static const double p_high = 0.97575;
48 
49  double x;
50  if ((p_low <= p) && (p <= p_high)) {
51  double q = p - 0.5;
52  double r = q * q;
53  x = (((((a[0] * r + a[1]) * r + a[2]) * r + a[3]) * r + a[4]) * r + a[5])
54  * q
55  / (((((b[0] * r + b[1]) * r + b[2]) * r + b[3]) * r + b[4]) * r + 1.0);
56  } else if (p < p_low) {
57  double q = std::sqrt(-2.0 * std::log(p));
58  x = (((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5])
59  / ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1.0);
60  } else {
61  double q = std::sqrt(-2.0 * log1m(p));
62  x = -(((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5])
63  / ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1.0);
64  }
65 
66  if (x < 37.6) { // gradient blows up past here
67  double e = Phi(x) - p;
68  double u = e * SQRT_2_TIMES_SQRT_PI * std::exp(0.5 * x * x);
69  x -= u / (1.0 + 0.5 * x * u);
70  }
71 
72  return x;
73 }
74 
75 } // namespace math
76 } // namespace stan
77 #endif
fvar< T > inv_Phi(const fvar< T > &p)
Definition: inv_Phi.hpp:14
void check_bounded(const char *function, const char *name, const T_y &y, const T_low &low, const T_high &high)
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
const double SQRT_2_TIMES_SQRT_PI
Definition: constants.hpp:136
const double a
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
Float_t d
Definition: plot.C:236
fvar< T > Phi(const fvar< T > &x)
Definition: Phi.hpp:12
double e()
Definition: constants.hpp:89
const hit & b
Definition: hits.cxx:21
TRandom3 r(0)
const double INFTY
Definition: constants.hpp:43
const double NEGATIVE_INFTY
Definition: constants.hpp:48
fvar< T > log1m(const fvar< T > &x)
Definition: log1m.hpp:12