Automatic Differentiation
 
Loading...
Searching...
No Matches
std_normal_log_qf.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_PROB_STD_NORMAL_LOG_QF_HPP
2#define STAN_MATH_PRIM_PROB_STD_NORMAL_LOG_QF_HPP
3
14#include <cmath>
15
16namespace stan {
17namespace math {
18
26inline double std_normal_log_qf(double log_p) {
27 check_not_nan("std_normal_log_qf", "Log probability variable", log_p);
28 check_less_or_equal("std_normal_log_qf", "Probability variable", log_p, 0);
29
30 if (log_p == NEGATIVE_INFTY) {
31 return NEGATIVE_INFTY;
32 }
33 if (log_p == 0) {
34 return INFTY;
35 }
36
37 static constexpr double log_a[8]
38 = {1.2199838032983212, 4.8914137334471356, 7.5865960847956080,
39 9.5274618535358388, 10.734698580862359, 11.116406781896242,
40 10.417226196842595, 7.8276718012189362};
41 static constexpr double log_b[8] = {0.,
42 3.7451021830139207,
43 6.5326064640478618,
44 8.5930788436817044,
45 9.9624069236663077,
46 10.579180688621286,
47 10.265665328832871,
48 8.5614962136628454};
49 static constexpr double log_c[8]
50 = {0.3530744474482423, 1.5326298343683388, 1.7525849400614634,
51 1.2941374937060454, 0.2393776640901312, -1.419724057885092,
52 -3.784340465764968, -7.163234779359426};
53 static constexpr double log_d[8] = {0.,
54 0.7193954734947205,
55 0.5166395879845317,
56 -0.371400933927844,
57 -1.909840708457214,
58 -4.186547581055928,
59 -7.509976771225415,
60 -20.67376157385924};
61 static constexpr double log_e[8]
62 = {1.8958048169567149, 1.6981417567726154, 0.5793212339927351,
63 -1.215503791936417, -3.629396584023968, -6.690500273261249,
64 -10.51540298415323, -15.41979457491781};
65 static constexpr double log_f[8] = {0.,
66 -0.511105318617135,
67 -1.988286302259815,
68 -4.208049039384857,
69 -7.147448611626374,
70 -10.89973190740069,
71 -15.76637472711685,
72 -33.82373901099482};
73
74 double log_q = log_p <= LOG_HALF ? log_diff_exp(LOG_HALF, log_p)
75 : log_diff_exp(log_p, LOG_HALF);
76 int log_q_sign = log_p <= LOG_HALF ? -1 : 1;
77 double log_r = log_q_sign == -1 ? log_p : log1m_exp(log_p);
78
79 if (stan::math::is_inf(log_r)) {
80 return 0;
81 }
82
83 double log_inner_r;
84 double log_pre_mult;
85 const double* num_ptr;
86 const double* den_ptr;
87
88 static constexpr double LOG_FIVE = LOG_TEN - LOG_TWO;
89 static constexpr double LOG_16 = LOG_TWO * 4;
90 static constexpr double LOG_425 = 6.0520891689244171729;
91 static constexpr double LOG_425_OVER_1000 = LOG_425 - LOG_TEN * 3;
92
93 if (log_q <= LOG_425_OVER_1000) {
94 log_inner_r = log_diff_exp(LOG_425_OVER_1000 * 2, log_q * 2);
95 log_pre_mult = log_q;
96 num_ptr = &log_a[0];
97 den_ptr = &log_b[0];
98 } else {
99 double log_temp_r = log(-log_r) / 2.0;
100 if (log_temp_r <= LOG_FIVE) {
101 log_inner_r = log_diff_exp(log_temp_r, LOG_16 - LOG_TEN);
102 num_ptr = &log_c[0];
103 den_ptr = &log_d[0];
104 } else {
105 log_inner_r = log_diff_exp(log_temp_r, LOG_FIVE);
106 num_ptr = &log_e[0];
107 den_ptr = &log_f[0];
108 }
109 log_pre_mult = 0.0;
110 }
111
112 // As computation requires evaluating r^8, this causes a loss of precision,
113 // even when on the log space. We can mitigate this by scaling the
114 // exponentiated result (dividing by 10), since the same scaling is applied
115 // to the numerator and denominator.
116 Eigen::VectorXd log_r_pow
117 = Eigen::ArrayXd::LinSpaced(8, 0, 7) * log_inner_r - LOG_TEN;
118 Eigen::Map<const Eigen::VectorXd> num_map(num_ptr, 8);
119 Eigen::Map<const Eigen::VectorXd> den_map(den_ptr, 8);
120 double log_result
121 = log_sum_exp(log_r_pow + num_map) - log_sum_exp(log_r_pow + den_map);
122 return log_q_sign * exp(log_pre_mult + log_result);
123}
124
134 template <typename T>
135 static inline auto fun(const T& x) {
136 return std_normal_log_qf(x);
137 }
138};
139
149template <
150 typename T,
153inline auto std_normal_log_qf(const T& x) {
155}
156
157} // namespace math
158} // namespace stan
159
160#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.
void check_less_or_equal(const char *function, const char *name, const T_y &y, const T_high &high, Idxs... idxs)
Throw an exception if y is not less than high.
static constexpr double LOG_HALF
The natural logarithm of 0.5, .
Definition constants.hpp:92
fvar< T > log1m_exp(const fvar< T > &x)
Return the natural logarithm of one minus the exponentiation of the specified argument.
Definition log1m_exp.hpp:23
static constexpr double LOG_TEN
The natural logarithm of 10, .
fvar< T > log(const fvar< T > &x)
Definition log.hpp:15
static constexpr double NEGATIVE_INFTY
Negative infinity.
Definition constants.hpp:51
static constexpr double LOG_TWO
The natural logarithm of 2, .
Definition constants.hpp:80
fvar< T > log_diff_exp(const fvar< T > &x1, const fvar< T > &x2)
void check_not_nan(const char *function, const char *name, const T_y &y)
Check if y is not NaN.
fvar< T > std_normal_log_qf(const fvar< T > &p)
int is_inf(const fvar< T > &x)
Returns 1 if the input's value is infinite and 0 otherwise.
Definition is_inf.hpp:21
static constexpr double INFTY
Positive infinity.
Definition constants.hpp:46
fvar< T > log_sum_exp(const fvar< T > &x1, const fvar< T > &x2)
fvar< T > exp(const fvar< T > &x)
Definition exp.hpp:13
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...
Structure to wrap std_normal_log_qf() so it can be vectorized.