1 #ifndef STAN_MATH_PRIM_PROB_EXP_MOD_NORMAL_CDF_HPP
2 #define STAN_MATH_PRIM_PROB_EXP_MOD_NORMAL_CDF_HPP
3
4 #include <stan/math/prim/meta.hpp>
5 #include <stan/math/prim/err.hpp>
6 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
7 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
8 #include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
9 #include <stan/math/prim/fun/constants.hpp>
10 #include <stan/math/prim/fun/erf.hpp>
11 #include <stan/math/prim/fun/exp.hpp>
12 #include <stan/math/prim/fun/inv.hpp>
13 #include <stan/math/prim/fun/is_inf.hpp>
14 #include <stan/math/prim/fun/max_size.hpp>
15 #include <stan/math/prim/fun/size.hpp>
16 #include <stan/math/prim/fun/size_zero.hpp>
17 #include <stan/math/prim/fun/square.hpp>
18 #include <stan/math/prim/fun/to_ref.hpp>
19 #include <stan/math/prim/fun/value_of.hpp>
20 #include <stan/math/prim/functor/operands_and_partials.hpp>
21 #include <cmath>
22
23 namespace stan {
24 namespace math {
25
26 template <typename T_y, typename T_loc, typename T_scale, typename T_inv_scale,
27 require_all_not_nonscalar_prim_or_rev_kernel_expression_t<
28 T_y, T_loc, T_scale, T_inv_scale>* = nullptr>
exp_mod_normal_cdf(const T_y & y,const T_loc & mu,const T_scale & sigma,const T_inv_scale & lambda)29 return_type_t<T_y, T_loc, T_scale, T_inv_scale> exp_mod_normal_cdf(
30 const T_y& y, const T_loc& mu, const T_scale& sigma,
31 const T_inv_scale& lambda) {
32 using T_partials_return = partials_return_t<T_y, T_loc, T_scale, T_inv_scale>;
33 using T_partials_array = Eigen::Array<T_partials_return, Eigen::Dynamic, 1>;
34 using T_y_ref = ref_type_if_t<!is_constant<T_y>::value, T_y>;
35 using T_mu_ref = ref_type_if_t<!is_constant<T_loc>::value, T_loc>;
36 using T_sigma_ref = ref_type_if_t<!is_constant<T_scale>::value, T_scale>;
37 using T_lambda_ref
38 = ref_type_if_t<!is_constant<T_inv_scale>::value, T_inv_scale>;
39 static const char* function = "exp_mod_normal_cdf";
40 check_consistent_sizes(function, "Random variable", y, "Location parameter",
41 mu, "Scale parameter", sigma, "Inv_scale paramter",
42 lambda);
43 T_y_ref y_ref = y;
44 T_mu_ref mu_ref = mu;
45 T_sigma_ref sigma_ref = sigma;
46 T_lambda_ref lambda_ref = lambda;
47
48 decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref));
49 decltype(auto) mu_val = to_ref(as_value_column_array_or_scalar(mu_ref));
50 decltype(auto) sigma_val = to_ref(as_value_column_array_or_scalar(sigma_ref));
51 decltype(auto) lambda_val
52 = to_ref(as_value_column_array_or_scalar(lambda_ref));
53
54 check_not_nan(function, "Random variable", y_val);
55 check_finite(function, "Location parameter", mu_val);
56 check_positive_finite(function, "Scale parameter", sigma_val);
57 check_positive_finite(function, "Inv_scale parameter", lambda_val);
58
59 if (size_zero(y, mu, sigma, lambda)) {
60 return 1.0;
61 }
62
63 operands_and_partials<T_y_ref, T_mu_ref, T_sigma_ref, T_lambda_ref>
64 ops_partials(y_ref, mu_ref, sigma_ref, lambda_ref);
65
66 using T_y_val_scalar = scalar_type_t<decltype(y_val)>;
67 if (is_vector<T_y>::value) {
68 if ((forward_as<Eigen::Array<T_y_val_scalar, Eigen::Dynamic, 1>>(y_val)
69 == NEGATIVE_INFTY)
70 .any()) {
71 return ops_partials.build(0.0);
72 }
73 } else {
74 if (forward_as<T_y_val_scalar>(y_val) == NEGATIVE_INFTY) {
75 return ops_partials.build(0.0);
76 }
77 }
78
79 const auto& inv_sigma
80 = to_ref_if<!is_constant_all<T_y, T_loc, T_scale>::value>(inv(sigma_val));
81 const auto& diff = to_ref(y_val - mu_val);
82 const auto& v = to_ref(lambda_val * sigma_val);
83 const auto& scaled_diff = to_ref(diff * INV_SQRT_TWO * inv_sigma);
84 const auto& scaled_diff_diff
85 = to_ref_if<!is_constant_all<T_y, T_loc, T_scale, T_inv_scale>::value>(
86 scaled_diff - v * INV_SQRT_TWO);
87 const auto& erf_calc = to_ref(0.5 * (1 + erf(scaled_diff_diff)));
88
89 const auto& exp_term
90 = to_ref_if<!is_constant_all<T_y, T_loc, T_scale, T_inv_scale>::value>(
91 exp(0.5 * square(v) - lambda_val * diff));
92 const auto& cdf_n
93 = to_ref(0.5 + 0.5 * erf(scaled_diff) - exp_term * erf_calc);
94
95 T_partials_return cdf(1.0);
96 if (is_vector<decltype(cdf_n)>::value) {
97 cdf = forward_as<T_partials_array>(cdf_n).prod();
98 } else {
99 cdf = forward_as<T_partials_return>(cdf_n);
100 }
101
102 if (!is_constant_all<T_y, T_loc, T_scale, T_inv_scale>::value) {
103 const auto& exp_term_2
104 = to_ref_if<(!is_constant_all<T_y, T_loc, T_scale>::value
105 && !is_constant_all<T_inv_scale>::value)>(
106 exp(-square(scaled_diff_diff)));
107 if (!is_constant_all<T_y, T_loc, T_scale>::value) {
108 constexpr bool need_deriv_refs = !is_constant_all<T_y, T_loc>::value
109 && !is_constant_all<T_scale>::value;
110 const auto& deriv_1
111 = to_ref_if<need_deriv_refs>(lambda_val * exp_term * erf_calc);
112 const auto& deriv_2 = to_ref_if<need_deriv_refs>(
113 INV_SQRT_TWO_PI * exp_term * exp_term_2 * inv_sigma);
114 const auto& sq_scaled_diff = square(scaled_diff);
115 const auto& exp_m_sq_scaled_diff = exp(-sq_scaled_diff);
116 const auto& deriv_3 = to_ref_if<need_deriv_refs>(
117 INV_SQRT_TWO_PI * exp_m_sq_scaled_diff * inv_sigma);
118 if (!is_constant_all<T_y, T_loc>::value) {
119 const auto& deriv = to_ref_if<(!is_constant_all<T_loc>::value
120 && !is_constant_all<T_y>::value)>(
121 cdf * (deriv_1 - deriv_2 + deriv_3) / cdf_n);
122 if (!is_constant_all<T_y>::value) {
123 ops_partials.edge1_.partials_ = deriv;
124 }
125 if (!is_constant_all<T_loc>::value) {
126 ops_partials.edge2_.partials_ = -deriv;
127 }
128 }
129 if (!is_constant_all<T_scale>::value) {
130 ops_partials.edge3_.partials_
131 = -cdf
132 * ((deriv_1 - deriv_2) * v
133 + (deriv_3 - deriv_2) * scaled_diff * SQRT_TWO)
134 / cdf_n;
135 }
136 }
137 if (!is_constant_all<T_inv_scale>::value) {
138 ops_partials.edge4_.partials_
139 = cdf * exp_term
140 * (INV_SQRT_TWO_PI * sigma_val * exp_term_2
141 - (v * sigma_val - diff) * erf_calc)
142 / cdf_n;
143 }
144 }
145 return ops_partials.build(cdf);
146 }
147
148 } // namespace math
149 } // namespace stan
150 #endif
151