| Stan Math Library 2.12.0 reverse mode automatic differentiation

- Main Page
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

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
4 #include <stan/math/prim/scal/fun/constants.hpp>
5 #include <stan/math/prim/scal/err/check_bounded.hpp>
6 #include <stan/math/prim/scal/fun/Phi.hpp>
7 #include <stan/math/prim/scal/fun/log1m.hpp>
8
9 namespace stan {
10 namespace math {
11
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,
36 -2.759285104469687e+02, 1.383577518672690e+02,
37 -3.066479806614716e+01, 2.506628277459239e+00
38 };
39 static const double b[5] = {
40 -5.447609879822406e+01, 1.615858368580409e+02,
41 -1.556989798598866e+02, 6.680131188771972e+01,
42 -1.328068155288572e+01
43 };
44 static const double c[6] = {
45 -7.784894002430293e-03, -3.223964580411365e-01,
46 -2.400758277161838e+00, -2.549732539343734e+00,
47 4.374664141464968e+00, 2.938163982698783e+00
48 };
49 static const double d[4] = {
50 7.784695709041462e-03, 3.224671290700398e-01,
51 2.445134137142996e+00, 3.754408661907416e+00
52 };
53
54 static const double p_low = 0.02425;
55 static const double p_high = 0.97575;
56
57 double x;
58 if ((p_low <= p) && (p <= p_high)) {
59 double q = p - 0.5;
60 double r = q * q;
61 x = (((((a[0]*r + a[1])*r + a[2])*r + a[3])*r + a[4])*r + a[5])*q
62 / (((((b[0]*r + b[1])*r + b[2])*r + b[3])*r + b[4])*r + 1.0);
63 } else if (p < p_low) {
64 double q = std::sqrt(-2.0*std::log(p));
65 x = (((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5])
66 / ((((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0);
67 } else {
68 double q = std::sqrt(-2.0 * log1m(p));
69 x = -(((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5])
70 / ((((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0);
71 }
72
73 if (x < 37.6) { // gradient blows up past here
74 double e = Phi(x) - p;
75 double u = e * SQRT_2_TIMES_SQRT_PI * std::exp(0.5 * x * x);
76 x -= u / (1.0 + 0.5 * x * u);
77 }
78
79 return x;
80 }
81
82 }
83 }
84 #endif
stan::math::sqrt
fvar< T > sqrt(const fvar< T > &x)
Definition: sqrt.hpp:14
stan::math::inv_Phi
fvar< T > inv_Phi(const fvar< T > &p)
Definition: inv_Phi.hpp:14
stan
Definition: log_sum_exp.hpp:8
stan::math::log
fvar< T > log(const fvar< T > &x)
Definition: log.hpp:14
stan::math::check_bounded
bool check_bounded(const char *function, const char *name, const T_y &y, const T_low &low, const T_high &high)
Return true if the value is between the low and high values, inclusively.
Definition: check_bounded.hpp:95
stan::math::SQRT_2_TIMES_SQRT_PI
const double SQRT_2_TIMES_SQRT_PI
Definition: constants.hpp:156
stan::math::exp
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
stan::math::Phi
fvar< T > Phi(const fvar< T > &x)
Definition: Phi.hpp:12
log1m.hpp
check_bounded.hpp
stan::math::e
double e()
Return the base of the natural logarithm.
Definition: constants.hpp:94
constants.hpp
stan::math::INFTY
const double INFTY
Positive infinity.
Definition: constants.hpp:43
stan::math::NEGATIVE_INFTY
const double NEGATIVE_INFTY
Negative infinity.
Definition: constants.hpp:49
Phi.hpp
stan::math::log1m
fvar< T > log1m(const fvar< T > &x)
Definition: log1m.hpp:15

[ Stan Home Page ] © 2011–2016, Stan Development Team.