1 #ifndef STAN_MATH_PRIM_PROB_NORMAL_LCDF_HPP
2 #define STAN_MATH_PRIM_PROB_NORMAL_LCDF_HPP
3
4 #include <stan/math/prim/meta.hpp>
5 #include <stan/math/prim/err.hpp>
6 #include <stan/math/prim/fun/constants.hpp>
7 #include <stan/math/prim/fun/erf.hpp>
8 #include <stan/math/prim/fun/erfc.hpp>
9 #include <stan/math/prim/fun/exp.hpp>
10 #include <stan/math/prim/fun/log.hpp>
11 #include <stan/math/prim/fun/log1p.hpp>
12 #include <stan/math/prim/fun/max_size.hpp>
13 #include <stan/math/prim/fun/scalar_seq_view.hpp>
14 #include <stan/math/prim/fun/size_zero.hpp>
15 #include <stan/math/prim/fun/square.hpp>
16 #include <stan/math/prim/fun/value_of.hpp>
17 #include <stan/math/prim/functor/operands_and_partials.hpp>
18 #include <cmath>
19 #include <limits>
20
21 namespace stan {
22 namespace math {
23
24 template <typename T_y, typename T_loc, typename T_scale,
25 require_all_not_nonscalar_prim_or_rev_kernel_expression_t<
26 T_y, T_loc, T_scale>* = nullptr>
normal_lcdf(const T_y & y,const T_loc & mu,const T_scale & sigma)27 inline return_type_t<T_y, T_loc, T_scale> normal_lcdf(const T_y& y,
28 const T_loc& mu,
29 const T_scale& sigma) {
30 using T_partials_return = partials_return_t<T_y, T_loc, T_scale>;
31 using std::exp;
32 using std::fabs;
33 using std::log;
34 using std::pow;
35 using std::sqrt;
36 using T_y_ref = ref_type_t<T_y>;
37 using T_mu_ref = ref_type_t<T_loc>;
38 using T_sigma_ref = ref_type_t<T_scale>;
39 static const char* function = "normal_lcdf";
40 check_consistent_sizes(function, "Random variable", y, "Location parameter",
41 mu, "Scale parameter", sigma);
42 T_y_ref y_ref = y;
43 T_mu_ref mu_ref = mu;
44 T_sigma_ref sigma_ref = sigma;
45 check_not_nan(function, "Random variable", y_ref);
46 check_finite(function, "Location parameter", mu_ref);
47 check_positive(function, "Scale parameter", sigma_ref);
48
49 if (size_zero(y, mu, sigma)) {
50 return 0;
51 }
52
53 T_partials_return cdf_log(0.0);
54 operands_and_partials<T_y_ref, T_mu_ref, T_sigma_ref> ops_partials(
55 y_ref, mu_ref, sigma_ref);
56
57 scalar_seq_view<T_y_ref> y_vec(y_ref);
58 scalar_seq_view<T_mu_ref> mu_vec(mu_ref);
59 scalar_seq_view<T_sigma_ref> sigma_vec(sigma_ref);
60 size_t N = max_size(y, mu, sigma);
61
62 for (size_t n = 0; n < N; n++) {
63 const T_partials_return y_dbl = y_vec.val(n);
64 const T_partials_return mu_dbl = mu_vec.val(n);
65 const T_partials_return sigma_dbl = sigma_vec.val(n);
66
67 const T_partials_return scaled_diff
68 = (y_dbl - mu_dbl) / (sigma_dbl * SQRT_TWO);
69
70 const T_partials_return sigma_sqrt2 = sigma_dbl * SQRT_TWO;
71 const T_partials_return x2 = square(scaled_diff);
72
73 // Rigorous numerical approximations are applied here to deal with values
74 // of |scaled_diff|>>0. This is needed to deal with rare base-rate
75 // logistic regression problems where it is useful to use an alternative
76 // link function instead.
77 //
78 // use erfc() instead of erf() in order to retain precision
79 // since for x>0 erfc()->0
80 if (scaled_diff > 0.0) {
81 // CDF(x) = 1/2 + 1/2erf(x) = 1 - 1/2erfc(x)
82 cdf_log += log1p(-0.5 * erfc(scaled_diff));
83 if (!is_not_nan(cdf_log)) {
84 cdf_log = 0;
85 }
86 } else if (scaled_diff > -20.0) {
87 // CDF(x) = 1/2 - 1/2erf(-x) = 1/2erfc(-x)
88 cdf_log += log(erfc(-scaled_diff)) + LOG_HALF;
89 } else if (10.0 * log(fabs(scaled_diff))
90 < log(std::numeric_limits<T_partials_return>::max())) {
91 // entering territory where erfc(-x)~0
92 // need to use direct numerical approximation of cdf_log instead
93 // the following based on W. J. Cody, Math. Comp. 23(107):631-638 (1969)
94 // CDF(x) = 1/2erfc(-x)
95 const T_partials_return x4 = pow(scaled_diff, 4);
96 const T_partials_return x6 = pow(scaled_diff, 6);
97 const T_partials_return x8 = pow(scaled_diff, 8);
98 const T_partials_return x10 = pow(scaled_diff, 10);
99 const T_partials_return temp_p
100 = 0.000658749161529837803157 + 0.0160837851487422766278 / x2
101 + 0.125781726111229246204 / x4 + 0.360344899949804439429 / x6
102 + 0.305326634961232344035 / x8 + 0.0163153871373020978498 / x10;
103 const T_partials_return temp_q
104 = -0.00233520497626869185443 - 0.0605183413124413191178 / x2
105 - 0.527905102951428412248 / x4 - 1.87295284992346047209 / x6
106 - 2.56852019228982242072 / x8 - 1.0 / x10;
107 cdf_log += LOG_HALF + log(INV_SQRT_PI + (temp_p / temp_q) / x2)
108 - log(-scaled_diff) - x2;
109 } else {
110 // scaled_diff^10 term will overflow
111 cdf_log = stan::math::negative_infinity();
112 }
113
114 if (!is_constant_all<T_y, T_loc, T_scale>::value) {
115 // compute partial derivatives
116 // based on analytic form given by:
117 // dln(CDF)/dx = exp(-x^2)/(sqrt(pi)*(1/2+erf(x)/2)
118 T_partials_return dncdf_log = 0.0;
119 T_partials_return t = 0.0;
120 T_partials_return t2 = 0.0;
121 T_partials_return t4 = 0.0;
122
123 // calculate using piecewise function
124 // (due to instability / inaccuracy in the various approximations)
125 if (scaled_diff > 2.9) {
126 // approximation derived from Abramowitz and Stegun (1964) 7.1.26
127 t = 1.0 / (1.0 + 0.3275911 * scaled_diff);
128 t2 = square(t);
129 t4 = pow(t, 4);
130 dncdf_log
131 = 1.0
132 / (SQRT_PI
133 * (exp(x2) - 0.254829592 + 0.284496736 * t - 1.421413741 * t2
134 + 1.453152027 * t2 * t - 1.061405429 * t4));
135 } else if (scaled_diff > 2.5) {
136 // in the trouble area where all of the standard numerical
137 // approximations are unstable - bridge the gap using Taylor
138 // expansions of the analytic function
139 // use Taylor expansion centred around x=2.7
140 t = scaled_diff - 2.7;
141 t2 = square(t);
142 t4 = pow(t, 4);
143 dncdf_log = 0.0003849882382 - 0.002079084702 * t + 0.005229340880 * t2
144 - 0.008029540137 * t2 * t + 0.008232190507 * t4
145 - 0.005692364250 * t4 * t + 0.002399496363 * pow(t, 6);
146 } else if (scaled_diff > 2.1) {
147 // use Taylor expansion centred around x=2.3
148 t = scaled_diff - 2.3;
149 t2 = square(t);
150 t4 = pow(t, 4);
151 dncdf_log = 0.002846135439 - 0.01310032351 * t + 0.02732189391 * t2
152 - 0.03326906904 * t2 * t + 0.02482478940 * t4
153 - 0.009883071924 * t4 * t - 0.0002771362254 * pow(t, 6);
154 } else if (scaled_diff > 1.5) {
155 // use Taylor expansion centred around x=1.85
156 t = scaled_diff - 1.85;
157 t2 = square(t);
158 t4 = pow(t, 4);
159 dncdf_log = 0.01849212058 - 0.06876280470 * t + 0.1099906382 * t2
160 - 0.09274533184 * t2 * t + 0.03543327418 * t4
161 + 0.005644855518 * t4 * t - 0.01111434424 * pow(t, 6);
162 } else if (scaled_diff > 0.8) {
163 // use Taylor expansion centred around x=1.15
164 t = scaled_diff - 1.15;
165 t2 = square(t);
166 t4 = pow(t, 4);
167 dncdf_log = 0.1585747034 - 0.3898677543 * t + 0.3515963775 * t2
168 - 0.09748053605 * t2 * t - 0.04347986191 * t4
169 + 0.02182506378 * t4 * t + 0.01074751427 * pow(t, 6);
170 } else if (scaled_diff > 0.1) {
171 // use Taylor expansion centred around x=0.45
172 t = scaled_diff - 0.45;
173 t2 = square(t);
174 t4 = pow(t, 4);
175 dncdf_log = 0.6245634904 - 0.9521866949 * t + 0.3986215682 * t2
176 + 0.04700850676 * t2 * t - 0.03478651979 * t4
177 - 0.01772675404 * t4 * t + 0.0006577254811 * pow(t, 6);
178 } else if (10.0 * log(fabs(scaled_diff))
179 < log(std::numeric_limits<T_partials_return>::max())) {
180 // approximation derived from Abramowitz and Stegun (1964) 7.1.26
181 // use fact that erf(x)=-erf(-x)
182 // Abramowitz and Stegun define this for -inf<x<0 but seems to be
183 // accurate for -inf<x<0.1
184 t = 1.0 / (1.0 - 0.3275911 * scaled_diff);
185 t2 = square(t);
186 t4 = pow(t, 4);
187 dncdf_log
188 = 2.0
189 / (SQRT_PI
190 * (0.254829592 * t - 0.284496736 * t2 + 1.421413741 * t2 * t
191 - 1.453152027 * t4 + 1.061405429 * t4 * t));
192 // check if we need to add a correction term
193 // (from cubic fit of residuals)
194 if (scaled_diff < -29.0) {
195 dncdf_log += 0.0015065154280332 * x2
196 - 0.3993154819705530 * scaled_diff - 4.2919418242931700;
197 } else if (scaled_diff < -17.0) {
198 dncdf_log += 0.0001263257217272 * x2 * scaled_diff
199 + 0.0123586859488623 * x2
200 - 0.0860505264736028 * scaled_diff - 1.252783383752970;
201 } else if (scaled_diff < -7.0) {
202 dncdf_log += 0.000471585349920831 * x2 * scaled_diff
203 + 0.0296839305424034 * x2
204 + 0.207402143352332 * scaled_diff + 0.425316974683324;
205 } else if (scaled_diff < -3.9) {
206 dncdf_log += -0.0006972280656443 * x2 * scaled_diff
207 + 0.0068218494628567 * x2
208 + 0.0585761964460277 * scaled_diff + 0.1034397670201370;
209 } else if (scaled_diff < -2.1) {
210 dncdf_log += -0.0018742199480885 * x2 * scaled_diff
211 - 0.0097119598291202 * x2
212 - 0.0170137970924080 * scaled_diff - 0.0100428567412041;
213 }
214 } else {
215 dncdf_log = stan::math::positive_infinity();
216 }
217
218 if (!is_constant_all<T_y>::value) {
219 ops_partials.edge1_.partials_[n] += dncdf_log / sigma_sqrt2;
220 }
221 if (!is_constant_all<T_loc>::value) {
222 ops_partials.edge2_.partials_[n] -= dncdf_log / sigma_sqrt2;
223 }
224 if (!is_constant_all<T_scale>::value) {
225 ops_partials.edge3_.partials_[n] -= dncdf_log * scaled_diff / sigma_dbl;
226 }
227 }
228 }
229 return ops_partials.build(cdf_log);
230 }
231
232 } // namespace math
233 } // namespace stan
234 #endif
235