Automatic Differentiation
 
Loading...
Searching...
No Matches
inv_Phi.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_FUN_INV_PHI_HPP
2#define STAN_MATH_PRIM_FUN_INV_PHI_HPP
3
11#include <cmath>
12
13namespace stan {
14namespace math {
15
16namespace internal {
23const int BIGINT = 2000000000;
24
31inline double inv_Phi_lambda(double p) {
32 check_bounded("inv_Phi", "Probability variable", p, 0, 1);
33
34 if (p < 8e-311) {
35 return NEGATIVE_INFTY;
36 }
37 if (p == 1) {
38 return INFTY;
39 }
40
41 static constexpr double a[8]
42 = {3.3871328727963666080e+00, 1.3314166789178437745e+02,
43 1.9715909503065514427e+03, 1.3731693765509461125e+04,
44 4.5921953931549871457e+04, 6.7265770927008700853e+04,
45 3.3430575583588128105e+04, 2.5090809287301226727e+03};
46 static constexpr double b[7]
47 = {4.2313330701600911252e+01, 6.8718700749205790830e+02,
48 5.3941960214247511077e+03, 2.1213794301586595867e+04,
49 3.9307895800092710610e+04, 2.8729085735721942674e+04,
50 5.2264952788528545610e+03};
51 static constexpr double c[8]
52 = {1.42343711074968357734e+00, 4.63033784615654529590e+00,
53 5.76949722146069140550e+00, 3.64784832476320460504e+00,
54 1.27045825245236838258e+00, 2.41780725177450611770e-01,
55 2.27238449892691845833e-02, 7.74545014278341407640e-04};
56 static constexpr double d[7]
57 = {2.05319162663775882187e+00, 1.67638483018380384940e+00,
58 6.89767334985100004550e-01, 1.48103976427480074590e-01,
59 1.51986665636164571966e-02, 5.47593808499534494600e-04,
60 1.05075007164441684324e-09};
61 static constexpr double e[8]
62 = {6.65790464350110377720e+00, 5.46378491116411436990e+00,
63 1.78482653991729133580e+00, 2.96560571828504891230e-01,
64 2.65321895265761230930e-02, 1.24266094738807843860e-03,
65 2.71155556874348757815e-05, 2.01033439929228813265e-07};
66 static constexpr double f[7]
67 = {5.99832206555887937690e-01, 1.36929880922735805310e-01,
68 1.48753612908506148525e-02, 7.86869131145613259100e-04,
69 1.84631831751005468180e-05, 1.42151175831644588870e-07,
70 2.04426310338993978564e-15};
71
72 double q = p - 0.5;
73 double r;
74 double val;
75
76 if (std::fabs(q) <= .425) {
77 r = .180625 - square(q);
78 return q
79 * (((((((a[7] * r + a[6]) * r + a[5]) * r + a[4]) * r + a[3]) * r
80 + a[2])
81 * r
82 + a[1])
83 * r
84 + a[0])
85 / (((((((b[6] * r + b[5]) * r + b[4]) * r + b[3]) * r + b[2]) * r
86 + b[1])
87 * r
88 + b[0])
89 * r
90 + 1.0);
91 } else {
92 r = q < 0 ? p : 1 - p;
93
94 if (r <= 0)
95 return 0;
96
97 r = std::sqrt(-std::log(r));
98
99 if (r <= 5.0) {
100 r += -1.6;
101 val = (((((((c[7] * r + c[6]) * r + c[5]) * r + c[4]) * r + c[3]) * r
102 + c[2])
103 * r
104 + c[1])
105 * r
106 + c[0])
107 / (((((((d[6] * r + d[5]) * r + d[4]) * r + d[3]) * r + d[2]) * r
108 + d[1])
109 * r
110 + d[0])
111 * r
112 + 1.0);
113 } else {
114 r -= 5.0;
115 val = (((((((e[7] * r + e[6]) * r + e[5]) * r + e[4]) * r + e[3]) * r
116 + e[2])
117 * r
118 + e[1])
119 * r
120 + e[0])
121 / (((((((f[6] * r + f[5]) * r + f[4]) * r + f[3]) * r + f[2]) * r
122 + f[1])
123 * r
124 + f[0])
125 * r
126 + 1.0);
127 }
128 if (q < 0.0)
129 return -val;
130 }
131 return val;
132}
133} // namespace internal
134
147inline double inv_Phi(double p) {
148 return p >= 0.9999 ? -internal::inv_Phi_lambda(
151}
152
162 template <typename T>
163 static inline auto fun(const T& x) {
164 return inv_Phi(x);
165 }
166};
167
176template <
177 typename T,
180inline auto inv_Phi(const T& x) {
182}
183
184} // namespace math
185} // namespace stan
186
187#endif
require_all_not_t< is_nonscalar_prim_or_rev_kernel_expression< std::decay_t< Types > >... > require_all_not_nonscalar_prim_or_rev_kernel_expression_t
Require none of the types satisfy is_nonscalar_prim_or_rev_kernel_expression.
require_not_t< is_var_matrix< std::decay_t< T > > > require_not_var_matrix_t
Require type does not satisfy is_var_matrix.
double inv_Phi_lambda(double p)
The inverse of the unit normal cumulative distribution function.
Definition inv_Phi.hpp:31
const int BIGINT
The largest integer that protects against floating point errors for the inv_Phi function.
Definition inv_Phi.hpp:23
void check_bounded(const char *function, const char *name, const T_y &y, const T_low &low, const T_high &high)
Check if the value is between the low and high values, inclusively.
static constexpr double e()
Return the base of the natural logarithm.
Definition constants.hpp:20
static constexpr double NEGATIVE_INFTY
Negative infinity.
Definition constants.hpp:51
fvar< T > inv_Phi(const fvar< T > &p)
Definition inv_Phi.hpp:15
static constexpr double INFTY
Positive infinity.
Definition constants.hpp:46
fvar< T > square(const fvar< T > &x)
Definition square.hpp:12
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...
Base template class for vectorization of unary scalar functions defined by a template class F to a sc...
static auto fun(const T &x)
Definition inv_Phi.hpp:163
Structure to wrap inv_Phi() so it can be vectorized.
Definition inv_Phi.hpp:161