1 #ifndef STAN_MATH_PRIM_PROB_LOGLOGISTIC_CDF_HPP
2 #define STAN_MATH_PRIM_PROB_LOGLOGISTIC_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/exp.hpp>
10 #include <stan/math/prim/fun/log.hpp>
11 #include <stan/math/prim/fun/max_size.hpp>
12 #include <stan/math/prim/fun/size.hpp>
13 #include <stan/math/prim/fun/size_zero.hpp>
14 #include <stan/math/prim/fun/to_ref.hpp>
15 #include <stan/math/prim/fun/value_of.hpp>
16 #include <stan/math/prim/functor/operands_and_partials.hpp>
17 #include <stan/math/prim/fun/promote_scalar.hpp>
18 #include <cmath>
19 #include <iostream>
20 
21 namespace stan {
22 namespace math {
23 
24 /** \ingroup prob_dists
25  * The loglogistic cumulative distribution function for the specified
26  * scalar(s) given the specified scales(s) and shape(s). y, alpha, or
27  * beta can each be either a scalar or a vector. Any vector inputs
28  * must be the same length.
29  *
30  *
31  * @tparam T_y type of scalar.
32  * @tparam T_scale type of scale parameter.
33  * @tparam T_shape type of shape parameter.
34  * @param y (Sequence of) scalar(s).
35  * @param alpha (Sequence of) scale parameter(s)
36  * for the loglogistic distribution.
37  * @param beta (Sequence of) shape parameter(s) for the
38  * loglogistic distribution.
39  * @return The loglogistic cdf evaluated at the specified arguments.
40  * @throw std::domain_error if any of the inputs are not positive or
41  * if and of the parameters are not finite.
42  */
43 template <typename T_y, typename T_scale, typename T_shape,
44           require_all_not_nonscalar_prim_or_rev_kernel_expression_t<
45               T_y, T_scale, T_shape>* = nullptr>
loglogistic_cdf(const T_y & y,const T_scale & alpha,const T_shape & beta)46 return_type_t<T_y, T_scale, T_shape> loglogistic_cdf(const T_y& y,
47                                                      const T_scale& alpha,
48                                                      const T_shape& beta) {
49   using T_partials_return = partials_return_t<T_y, T_scale, T_shape>;
50   using T_y_ref = ref_type_t<T_y>;
51   using T_alpha_ref = ref_type_t<T_scale>;
52   using T_beta_ref = ref_type_t<T_shape>;
53   using std::pow;
54   static const char* function = "loglogistic_cdf";
55   check_consistent_sizes(function, "Random variable", y, "Scale parameter",
56                          alpha, "Shape parameter", beta);
57   T_y_ref y_ref = y;
58   T_alpha_ref alpha_ref = alpha;
59   T_beta_ref beta_ref = beta;
60 
61   decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref));
62   decltype(auto) alpha_val = to_ref(as_value_column_array_or_scalar(alpha_ref));
63   decltype(auto) beta_val = to_ref(as_value_column_array_or_scalar(beta_ref));
64 
65   check_nonnegative(function, "Random variable", y_val);
66   check_positive_finite(function, "Scale parameter", alpha_val);
67   check_positive_finite(function, "Shape parameter", beta_val);
68 
69   if (size_zero(y, alpha, beta)) {
70     return 1.0;
71   }
72 
73   operands_and_partials<T_y_ref, T_alpha_ref, T_beta_ref> ops_partials(
74       y_ref, alpha_ref, beta_ref);
75 
76   if (sum(promote_scalar<int>(y_val == 0))) {
77     return ops_partials.build(0.0);
78   }
79 
80   const auto& alpha_div_y
81       = to_ref_if<!is_constant_all<T_shape>::value>(alpha_val / y_val);
82   const auto& alpha_div_y_pow_beta
83       = to_ref_if<!is_constant_all<T_y, T_scale, T_shape>::value>(
84           pow(alpha_div_y, beta_val));
85   const auto& prod_all
86       = to_ref_if<!is_constant_all<T_y, T_scale, T_shape>::value>(
87           1 / (1 + alpha_div_y_pow_beta));
88 
89   T_partials_return cdf = prod(prod_all);
90 
91   if (!is_constant_all<T_y, T_scale, T_shape>::value) {
92     const auto& prod_all_sq = to_ref_if<!is_constant_all<T_y>::value
93                                             + !is_constant_all<T_scale>::value
94                                             + !is_constant_all<T_shape>::value
95                                         >= 2>(square(prod_all));
96     const auto& cdf_div_elt = to_ref_if<!is_constant_all<T_y>::value
97                                             + !is_constant_all<T_scale>::value
98                                             + !is_constant_all<T_shape>::value
99                                         >= 2>(cdf / prod_all);
100     if (!is_constant_all<T_y, T_scale>::value) {
101       const auto& alpha_div_times_beta = to_ref_if<
102           !is_constant_all<T_y>::value + !is_constant_all<T_scale>::value == 2>(
103           alpha_div_y_pow_beta * beta_val);
104       if (!is_constant_all<T_y>::value) {
105         const auto& y_deriv = alpha_div_times_beta / y_val * prod_all_sq;
106         ops_partials.edge1_.partials_ = y_deriv * cdf_div_elt;
107       }
108       if (!is_constant_all<T_scale>::value) {
109         const auto& alpha_deriv
110             = -alpha_div_times_beta / alpha_val * prod_all_sq;
111         ops_partials.edge2_.partials_ = alpha_deriv * cdf_div_elt;
112       }
113     }
114     if (!is_constant_all<T_shape>::value) {
115       const auto& beta_deriv
116           = -multiply_log(alpha_div_y_pow_beta, alpha_div_y) * prod_all_sq;
117       ops_partials.edge3_.partials_ = beta_deriv * cdf_div_elt;
118     }
119   }
120 
121   return ops_partials.build(cdf);
122 }
123 
124 }  // namespace math
125 }  // namespace stan
126 #endif
127