1 #ifndef STAN_MATH_REV_FUN_PHI_APPROX_HPP
2 #define STAN_MATH_REV_FUN_PHI_APPROX_HPP
3
4 #include <stan/math/rev/meta.hpp>
5 #include <stan/math/rev/core.hpp>
6 #include <stan/math/prim/fun/inv_logit.hpp>
7
8 namespace stan {
9 namespace math {
10
11 /**
12 * Approximation of the unit normal CDF for variables (stan).
13 *
14 * http://www.jiem.org/index.php/jiem/article/download/60/27
15 *
16 *
17 \f[
18 \mbox{Phi\_approx}(x) =
19 \begin{cases}
20 \Phi_{\mbox{\footnotesize approx}}(x) & \mbox{if } -\infty\leq x\leq \infty
21 \\[6pt] \textrm{NaN} & \mbox{if } x = \textrm{NaN} \end{cases} \f]
22
23 \f[
24 \frac{\partial\, \mbox{Phi\_approx}(x)}{\partial x} =
25 \begin{cases}
26 \frac{\partial\, \Phi_{\mbox{\footnotesize approx}}(x)}{\partial x}
27 & \mbox{if } -\infty\leq x\leq \infty \\[6pt]
28 \textrm{NaN} & \mbox{if } x = \textrm{NaN}
29 \end{cases}
30 \f]
31
32 \f[
33 \Phi_{\mbox{\footnotesize approx}}(x) = \mbox{logit}^{-1}(0.07056 \,
34 x^3 + 1.5976 \, x)
35 \f]
36
37 \f[
38 \frac{\partial \, \Phi_{\mbox{\footnotesize approx}}(x)}{\partial x}
39 = -\Phi_{\mbox{\footnotesize approx}}^2(x)
40 e^{-0.07056x^3 - 1.5976x}(-0.21168x^2-1.5976)
41 \f]
42 *
43 * @param a Variable argument.
44 * @return The corresponding unit normal cdf approximation.
45 */
Phi_approx(const var & a)46 inline var Phi_approx(const var& a) {
47 double av_squared = a.val() * a.val();
48 double f = inv_logit(0.07056 * a.val() * av_squared + 1.5976 * a.val());
49 double da = f * (1 - f) * (3.0 * 0.07056 * av_squared + 1.5976);
50 return make_callback_var(
51 f, [a, da](auto& vi) mutable { a.adj() += vi.adj() * da; });
52 }
53
54 template <typename T, require_var_matrix_t<T>* = nullptr>
Phi_approx(const T & a)55 inline auto Phi_approx(const T& a) {
56 arena_t<value_type_t<T>> f(a.rows(), a.cols());
57 arena_t<value_type_t<T>> da(a.rows(), a.cols());
58 for (Eigen::Index j = 0; j < a.cols(); ++j) {
59 for (Eigen::Index i = 0; i < a.rows(); ++i) {
60 const auto a_val = a.val().coeff(i, j);
61 const auto av_squared = a_val * a_val;
62 f.coeffRef(i, j) = inv_logit(0.07056 * a_val * av_squared
63 + 1.5976 * a.val().coeff(i, j));
64 da.coeffRef(i, j) = f.coeff(i, j) * (1 - f.coeff(i, j))
65 * (3.0 * 0.07056 * av_squared + 1.5976);
66 }
67 }
68 return make_callback_var(f, [a, da](auto& vi) mutable {
69 a.adj().array() += vi.adj().array() * da.array();
70 });
71 }
72
73 } // namespace math
74 } // namespace stan
75 #endif
76