Automatic Differentiation
 
Loading...
Searching...
No Matches
normal_cdf.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_PROB_NORMAL_CDF_HPP
2#define STAN_MATH_PRIM_PROB_NORMAL_CDF_HPP
3
16#include <cmath>
17
18namespace stan {
19namespace math {
20
35template <typename T_y, typename T_loc, typename T_scale,
37 T_y, T_loc, T_scale>* = nullptr>
39 const T_loc& mu,
40 const T_scale& sigma) {
41 using T_partials_return = partials_return_t<T_y, T_loc, T_scale>;
42 using std::exp;
43 using T_y_ref = ref_type_t<T_y>;
44 using T_mu_ref = ref_type_t<T_loc>;
45 using T_sigma_ref = ref_type_t<T_scale>;
46 static constexpr const char* function = "normal_cdf";
47 check_consistent_sizes(function, "Random variable", y, "Location parameter",
48 mu, "Scale parameter", sigma);
49 T_y_ref y_ref = y;
50 T_mu_ref mu_ref = mu;
51 T_sigma_ref sigma_ref = sigma;
52 check_not_nan(function, "Random variable", y_ref);
53 check_finite(function, "Location parameter", mu_ref);
54 check_positive(function, "Scale parameter", sigma_ref);
55
56 if (size_zero(y, mu, sigma)) {
57 return 1.0;
58 }
59
60 T_partials_return cdf(1.0);
61 auto ops_partials = make_partials_propagator(y_ref, mu_ref, sigma_ref);
62
63 scalar_seq_view<T_y_ref> y_vec(y_ref);
64 scalar_seq_view<T_mu_ref> mu_vec(mu_ref);
65 scalar_seq_view<T_sigma_ref> sigma_vec(sigma_ref);
66 size_t N = max_size(y, mu, sigma);
67
68 for (size_t n = 0; n < N; n++) {
69 const T_partials_return y_dbl = y_vec.val(n);
70 const T_partials_return mu_dbl = mu_vec.val(n);
71 const T_partials_return sigma_dbl = sigma_vec.val(n);
72 const T_partials_return scaled_diff
73 = (y_dbl - mu_dbl) / (sigma_dbl * SQRT_TWO);
74 T_partials_return cdf_n;
75 if (scaled_diff < -37.5 * INV_SQRT_TWO) {
76 cdf_n = 0.0;
77 } else if (scaled_diff < -5.0 * INV_SQRT_TWO) {
78 cdf_n = 0.5 * erfc(-scaled_diff);
79 } else if (scaled_diff > 8.25 * INV_SQRT_TWO) {
80 cdf_n = 1;
81 } else {
82 cdf_n = 0.5 * (1.0 + erf(scaled_diff));
83 }
84
85 cdf *= cdf_n;
86
88 const T_partials_return rep_deriv
89 = (scaled_diff < -37.5 * INV_SQRT_TWO)
90 ? 0.0
91 : INV_SQRT_TWO_PI * exp(-scaled_diff * scaled_diff)
92 / (cdf_n * sigma_dbl);
94 partials<0>(ops_partials)[n] += rep_deriv;
95 }
97 partials<1>(ops_partials)[n] -= rep_deriv;
98 }
100 partials<2>(ops_partials)[n] -= rep_deriv * scaled_diff * SQRT_TWO;
101 }
102 }
103 }
104
106 for (size_t n = 0; n < stan::math::size(y); ++n) {
107 partials<0>(ops_partials)[n] *= cdf;
108 }
109 }
111 for (size_t n = 0; n < stan::math::size(mu); ++n) {
112 partials<1>(ops_partials)[n] *= cdf;
113 }
114 }
116 for (size_t n = 0; n < stan::math::size(sigma); ++n) {
117 partials<2>(ops_partials)[n] *= cdf;
118 }
119 }
120 return ops_partials.build(cdf);
121}
122
123} // namespace math
124} // namespace stan
125#endif
scalar_seq_view provides a uniform sequence-like wrapper around either a scalar or a sequence of scal...
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.
return_type_t< T_y_cl, T_loc_cl, T_scale_cl > normal_cdf(const T_y_cl &y, const T_loc_cl &mu, const T_scale_cl &sigma)
Returns the normal cumulative distribution function for the given location, and scale.
typename return_type< Ts... >::type return_type_t
Convenience type for the return type of the specified template parameters.
int64_t size(const T &m)
Returns the size (number of the elements) of a matrix_cl or var_value<matrix_cl<T>>.
Definition size.hpp:19
bool size_zero(const T &x)
Returns 1 if input is of length 0, returns 0 otherwise.
Definition size_zero.hpp:19
fvar< T > erf(const fvar< T > &x)
Definition erf.hpp:15
static constexpr double INV_SQRT_TWO
The value of 1 over the square root of 2, .
static constexpr double INV_SQRT_TWO_PI
The value of 1 over the square root of , .
static constexpr double SQRT_TWO
The value of the square root of 2, .
void check_consistent_sizes(const char *)
Trivial no input case, this function is a no-op.
fvar< T > erfc(const fvar< T > &x)
Definition erfc.hpp:15
void check_finite(const char *function, const char *name, const T_y &y)
Return true if all values in y are finite.
void check_not_nan(const char *function, const char *name, const T_y &y)
Check if y is not NaN.
void check_positive(const char *function, const char *name, const T_y &y)
Check if y is positive.
int64_t max_size(const T1 &x1, const Ts &... xs)
Calculate the size of the largest input.
Definition max_size.hpp:20
auto make_partials_propagator(Ops &&... ops)
Construct an partials_propagator.
fvar< T > exp(const fvar< T > &x)
Definition exp.hpp:13
typename ref_type_if< true, T >::type ref_type_t
Definition ref_type.hpp:55
typename partials_return_type< Args... >::type partials_return_t
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...
Extends std::true_type when instantiated with zero or more template parameters, all of which extend t...