Skip to content

Commit 29843fe

Browse files
rok-cesnovarrok-cesnovar
authored
Merge pull request #2488 from gregorp90/feature/loglogistic_cdf
Added loglogistic cdf and rng
2 parents 5534e22 + b09414e commit 29843fe

5 files changed

Lines changed: 354 additions & 0 deletions

File tree

‎stan/math/prim/prob.hpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -175,7 +175,9 @@
175175
#include <stan/math/prim/prob/logistic_log.hpp>
176176
#include <stan/math/prim/prob/logistic_lpdf.hpp>
177177
#include <stan/math/prim/prob/logistic_rng.hpp>
178+
#include <stan/math/prim/prob/loglogistic_cdf.hpp>
178179
#include <stan/math/prim/prob/loglogistic_lpdf.hpp>
180+
#include <stan/math/prim/prob/loglogistic_rng.hpp>
179181
#include <stan/math/prim/prob/lognormal_ccdf_log.hpp>
180182
#include <stan/math/prim/prob/lognormal_cdf.hpp>
181183
#include <stan/math/prim/prob/lognormal_cdf_log.hpp>

‎stan/math/prim/prob/loglogistic_cdf.hpp‎

Lines changed: 126 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,126 @@
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>
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

‎stan/math/prim/prob/loglogistic_rng.hpp‎

Lines changed: 63 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,63 @@
1+
#ifndef STAN_MATH_PRIM_PROB_LOGLOGISTIC_RNG_HPP
2+
#define STAN_MATH_PRIM_PROB_LOGLOGISTIC_RNG_HPP
3+
4+
#include <stan/math/prim/meta.hpp>
5+
#include <stan/math/prim/err.hpp>
6+
#include <stan/math/prim/fun/scalar_seq_view.hpp>
7+
#include <stan/math/prim/fun/max_size.hpp>
8+
#include <boost/random/uniform_01.hpp>
9+
#include <boost/random/variate_generator.hpp>
10+
11+
namespace stan {
12+
namespace math {
13+
14+
/** \ingroup prob_dists
15+
* Return a loglogistic random variate for the given scale and
16+
* shape parameters using the specified random number generator.
17+
*
18+
* alpha and beta can each be a scalar or a one-dimensional container. Any
19+
* non-scalar inputs must be the same size.
20+
*
21+
* @tparam T_scale type of scale parameter
22+
* @tparam T_shape type of shape parameter
23+
* @tparam RNG type of random number generator
24+
* @param alpha (Sequence of) positive scale parameter(s)
25+
* @param beta (Sequence of) positive shape parameter(s)
26+
* @param rng random number generator
27+
* @return (Sequence of) loglogistic random variate(s)
28+
* @throw std::domain_error if alpha or beta are nonpositive
29+
* @throw std::invalid_argument if non-scalar arguments are of different
30+
* sizes
31+
*/
32+
template <typename T_scale, typename T_shape, class RNG>
33+
inline typename VectorBuilder<true, double, T_scale, T_shape>::type
34+
loglogistic_rng(const T_scale& alpha, const T_shape& beta, RNG& rng) {
35+
using boost::uniform_01;
36+
using boost::variate_generator;
37+
using std::pow;
38+
using T_alpha_ref = ref_type_t<T_scale>;
39+
using T_beta_ref = ref_type_t<T_shape>;
40+
static const char* function = "loglogistic_rng";
41+
check_consistent_sizes(function, "Scale parameter", alpha, "Shape Parameter",
42+
beta);
43+
T_alpha_ref alpha_ref = alpha;
44+
T_beta_ref beta_ref = beta;
45+
check_positive_finite(function, "Scale parameter", alpha_ref);
46+
check_positive_finite(function, "Shape parameter", beta_ref);
47+
48+
scalar_seq_view<T_alpha_ref> alpha_vec(alpha_ref);
49+
scalar_seq_view<T_beta_ref> beta_vec(beta_ref);
50+
size_t N = max_size(alpha, beta);
51+
VectorBuilder<true, double, T_scale, T_shape> output(N);
52+
53+
for (size_t n = 0; n < N; ++n) {
54+
variate_generator<RNG&, uniform_01<> > uniform01_rng(rng, uniform_01<>());
55+
const double tmp = uniform01_rng();
56+
output[n] = alpha_vec[n] * pow(tmp / (1 - tmp), 1 / beta_vec[n]);
57+
}
58+
return output.data();
59+
}
60+
61+
} // namespace math
62+
} // namespace stan
63+
#endif

‎test/prob/loglogistic/loglogistic_cdf_test.hpp‎

Lines changed: 86 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,86 @@
1+
// Arguments: Doubles, Doubles, Doubles
2+
#include <stan/math/prim.hpp>
3+
4+
using stan::math::var;
5+
using std::numeric_limits;
6+
using std::pow;
7+
using std::vector;
8+
9+
class AgradCdfLogistic : public AgradCdfTest {
10+
public:
11+
void valid_values(vector<vector<double> >& parameters, vector<double>& cdf) {
12+
vector<double> param(3);
13+
14+
param[0] = 2.0; // y
15+
param[1] = 3.0; // Scale
16+
param[2] = 2.0; // Shape
17+
parameters.push_back(param);
18+
cdf.push_back(0.30769230769230765388); // expected cdf
19+
20+
param[0] = 20.0; // y
21+
param[1] = 11.5; // Scale
22+
param[2] = 15.0; // Shape
23+
parameters.push_back(param);
24+
cdf.push_back(0.99975173823523311167); // expected cdf
25+
26+
param[0] = 0.0001; // y
27+
param[1] = 1.2; // Scale
28+
param[2] = 1.0; // Shape
29+
parameters.push_back(param);
30+
cdf.push_back(0.00008332638946754438); // expected cdf
31+
}
32+
33+
void invalid_values(vector<size_t>& index, vector<double>& value) {
34+
// y
35+
index.push_back(0U);
36+
value.push_back(-1.0);
37+
38+
index.push_back(0U);
39+
value.push_back(-numeric_limits<double>::infinity());
40+
41+
// alpha
42+
index.push_back(1U);
43+
value.push_back(-1);
44+
45+
index.push_back(1U);
46+
value.push_back(0);
47+
48+
index.push_back(1U);
49+
value.push_back(numeric_limits<double>::infinity());
50+
51+
// beta
52+
index.push_back(2U);
53+
value.push_back(-1.0);
54+
55+
index.push_back(2U);
56+
value.push_back(0.0);
57+
58+
index.push_back(2U);
59+
value.push_back(numeric_limits<double>::infinity());
60+
}
61+
62+
bool has_lower_bound() { return true; }
63+
64+
double lower_bound() { return 0.0; }
65+
66+
bool has_upper_bound() { return false; }
67+
68+
template <typename T_y, typename T_scale, typename T_shape, typename T3,
69+
typename T4, typename T5>
70+
stan::return_type_t<T_y, T_scale, T_shape> cdf(const T_y& y,
71+
const T_scale& alpha,
72+
const T_shape& beta, const T3&,
73+
const T4&, const T5&) {
74+
return stan::math::loglogistic_cdf(y, alpha, beta);
75+
}
76+
77+
template <typename T_y, typename T_scale, typename T_shape, typename T3,
78+
typename T4, typename T5>
79+
stan::return_type_t<T_y, T_scale, T_shape> cdf_function(const T_y& y,
80+
const T_scale& alpha,
81+
const T_shape& beta,
82+
const T3&, const T4&,
83+
const T5&) {
84+
return 1.0 / (1.0 + pow(alpha / y, beta));
85+
}
86+
};

‎test/unit/math/prim/prob/loglogistic_test.cpp‎

Lines changed: 77 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,77 @@
1+
#include <stan/math/prim.hpp>
2+
#include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
3+
#include <test/unit/math/prim/prob/util.hpp>
4+
#include <gtest/gtest.h>
5+
#include <boost/math/distributions.hpp>
6+
#include <boost/random/mersenne_twister.hpp>
7+
#include <limits>
8+
#include <vector>
9+
10+
class LoglogisticTestRig : public VectorRNGTestRig {
11+
public:
12+
LoglogisticTestRig()
13+
: VectorRNGTestRig(10000, 10, {2.5, 1.7, 0.2, 0.1, 2.0},
14+
{3, 2, 1, 5, 10, 6}, {-2.5, -1.7, -0.2, -0.1, 0.0},
15+
{-3, -2, -1, -4, -10, 0}, {0.1, 1.0, 2.5, 4.0},
16+
{1, 2, 3, 4}, {-2.7, -1.5, -0.5, 0.0},
17+
{-3, -2, -1, 0}) {}
18+
19+
template <typename T1, typename T2, typename T3, typename T_rng>
20+
auto generate_samples(const T1& alpha, const T2& beta, const T3& unused,
21+
T_rng& rng) const {
22+
return stan::math::loglogistic_rng(alpha, beta, rng);
23+
}
24+
};
25+
26+
double icdf(double x, double alpha, double beta) {
27+
return alpha * pow(x / (1 - x), 1 / beta);
28+
}
29+
30+
TEST(ProbDistributionsLoglogistic, errorCheck) {
31+
check_dist_throws_all_types(LoglogisticTestRig());
32+
}
33+
34+
TEST(ProbDistributionsLoglogistic, error_check) {
35+
boost::random::mt19937 rng;
36+
EXPECT_NO_THROW(stan::math::loglogistic_rng(10.0, 2.0, rng));
37+
38+
EXPECT_THROW(stan::math::loglogistic_rng(2.0, -1.0, rng), std::domain_error);
39+
EXPECT_THROW(stan::math::loglogistic_rng(-2.0, 1.0, rng), std::domain_error);
40+
EXPECT_THROW(
41+
stan::math::loglogistic_rng(10, stan::math::positive_infinity(), rng),
42+
std::domain_error);
43+
EXPECT_THROW(
44+
stan::math::loglogistic_rng(stan::math::positive_infinity(), 2, rng),
45+
std::domain_error);
46+
}
47+
48+
TEST(ProbDistributionsLoglogistic, test_sampling_icdf) {
49+
for (double p : {0.0, 0.1, 0.2, 0.5, 0.7, 0.9, 0.99}) {
50+
for (double alpha : {1.11, 0.13, 1.2, 4.67}) {
51+
for (double beta : {0.11, 1.33, 2.0, 3.2}) {
52+
double x = icdf(p, alpha, beta);
53+
EXPECT_FLOAT_EQ(stan::math::loglogistic_cdf(x, alpha, beta), p);
54+
}
55+
}
56+
}
57+
}
58+
59+
TEST(ProbDistributionsLoglogistic, chiSquareGoodnessFitTest) {
60+
boost::random::mt19937 rng;
61+
int N = 10000;
62+
int K = stan::math::round(2 * std::pow(N, 0.4));
63+
64+
std::vector<double> samples;
65+
for (int i = 0; i < N; ++i) {
66+
samples.push_back(stan::math::loglogistic_rng(1.2, 2.0, rng));
67+
}
68+
std::vector<double> quantiles;
69+
for (int i = 1; i < K; ++i) {
70+
double frac = static_cast<double>(i) / K;
71+
quantiles.push_back(icdf(frac, 1.2, 2.0));
72+
}
73+
quantiles.push_back(std::numeric_limits<double>::max());
74+
75+
// Assert that they match
76+
assert_matches_quantiles(samples, quantiles, 1e-6);
77+
}

0 commit comments

Comments
 (0)