Skip to content

Commit 18dbb5b

Browse files
bbbales2bbbales2
authored
Merge pull request #2136 from bstatcomp/generalize_weibull
generalize Weibull
2 parents 5aa6d2f + 3473789 commit 18dbb5b

File tree

6 files changed

+193
-591
lines changed

6 files changed

+193
-591
lines changed

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

+46-42
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88
#include <stan/math/prim/fun/max_size.hpp>
99
#include <stan/math/prim/fun/size.hpp>
1010
#include <stan/math/prim/fun/size_zero.hpp>
11+
#include <stan/math/prim/fun/to_ref.hpp>
1112
#include <stan/math/prim/fun/value_of.hpp>
1213
#include <stan/math/prim/functor/operands_and_partials.hpp>
1314
#include <cmath>
@@ -34,60 +35,63 @@ return_type_t<T_y, T_shape, T_scale> weibull_cdf(const T_y& y,
3435
const T_shape& alpha,
3536
const T_scale& sigma) {
3637
using T_partials_return = partials_return_t<T_y, T_shape, T_scale>;
37-
using std::exp;
38-
using std::log;
38+
using T_y_ref = ref_type_if_t<!is_constant<T_y>::value, T_y>;
39+
using T_alpha_ref = ref_type_if_t<!is_constant<T_shape>::value, T_shape>;
40+
using T_sigma_ref = ref_type_if_t<!is_constant<T_scale>::value, T_scale>;
3941
using std::pow;
4042
static const char* function = "weibull_cdf";
41-
check_nonnegative(function, "Random variable", y);
42-
check_positive_finite(function, "Shape parameter", alpha);
43-
check_positive_finite(function, "Scale parameter", sigma);
43+
44+
T_y_ref y_ref = y;
45+
T_alpha_ref alpha_ref = alpha;
46+
T_sigma_ref sigma_ref = sigma;
47+
48+
const auto& y_col = as_column_vector_or_scalar(y_ref);
49+
const auto& alpha_col = as_column_vector_or_scalar(alpha_ref);
50+
const auto& sigma_col = as_column_vector_or_scalar(sigma_ref);
51+
52+
const auto& y_arr = as_array_or_scalar(y_col);
53+
const auto& alpha_arr = as_array_or_scalar(alpha_col);
54+
const auto& sigma_arr = as_array_or_scalar(sigma_col);
55+
56+
ref_type_t<decltype(value_of(y_arr))> y_val = value_of(y_arr);
57+
ref_type_t<decltype(value_of(alpha_arr))> alpha_val = value_of(alpha_arr);
58+
ref_type_t<decltype(value_of(sigma_arr))> sigma_val = value_of(sigma_arr);
59+
60+
check_nonnegative(function, "Random variable", y_val);
61+
check_positive_finite(function, "Shape parameter", alpha_val);
62+
check_positive_finite(function, "Scale parameter", sigma_val);
4463

4564
if (size_zero(y, alpha, sigma)) {
4665
return 1.0;
4766
}
4867

49-
T_partials_return cdf(1.0);
50-
operands_and_partials<T_y, T_shape, T_scale> ops_partials(y, alpha, sigma);
68+
operands_and_partials<T_y_ref, T_alpha_ref, T_sigma_ref> ops_partials(
69+
y_ref, alpha_ref, sigma_ref);
5170

52-
scalar_seq_view<T_y> y_vec(y);
53-
scalar_seq_view<T_scale> sigma_vec(sigma);
54-
scalar_seq_view<T_shape> alpha_vec(alpha);
55-
size_t N = max_size(y, sigma, alpha);
56-
for (size_t n = 0; n < N; n++) {
57-
const T_partials_return y_dbl = value_of(y_vec[n]);
58-
const T_partials_return sigma_dbl = value_of(sigma_vec[n]);
59-
const T_partials_return alpha_dbl = value_of(alpha_vec[n]);
60-
const T_partials_return pow_n = pow(y_dbl / sigma_dbl, alpha_dbl);
61-
const T_partials_return exp_n = exp(-pow_n);
62-
const T_partials_return cdf_n = 1.0 - exp_n;
71+
constexpr bool any_derivs = !is_constant_all<T_y, T_shape, T_scale>::value;
72+
const auto& pow_n = to_ref_if<any_derivs>(pow(y_val / sigma_val, alpha_val));
73+
const auto& exp_n = to_ref_if<any_derivs>(exp(-pow_n));
74+
const auto& cdf_n = to_ref_if<any_derivs>(1 - exp_n);
6375

64-
cdf *= cdf_n;
76+
T_partials_return cdf = prod(cdf_n);
6577

66-
const T_partials_return rep_deriv = exp_n * pow_n / cdf_n;
67-
if (!is_constant_all<T_y>::value) {
68-
ops_partials.edge1_.partials_[n] += rep_deriv * alpha_dbl / y_dbl;
78+
if (any_derivs) {
79+
const auto& rep_deriv = to_ref_if<(!is_constant_all<T_y, T_scale>::value
80+
&& !is_constant_all<T_shape>::value)>(
81+
exp_n * pow_n * cdf / cdf_n);
82+
if (!is_constant_all<T_y, T_scale>::value) {
83+
const auto& deriv_y_sigma = to_ref_if<(
84+
!is_constant_all<T_y>::value && !is_constant_all<T_scale>::value)>(
85+
rep_deriv * alpha_val);
86+
if (!is_constant_all<T_y>::value) {
87+
ops_partials.edge1_.partials_ = deriv_y_sigma / y_val;
88+
}
89+
if (!is_constant_all<T_scale>::value) {
90+
ops_partials.edge3_.partials_ = -deriv_y_sigma / sigma_val;
91+
}
6992
}
7093
if (!is_constant_all<T_shape>::value) {
71-
ops_partials.edge2_.partials_[n] += rep_deriv * log(y_dbl / sigma_dbl);
72-
}
73-
if (!is_constant_all<T_scale>::value) {
74-
ops_partials.edge3_.partials_[n] -= rep_deriv * alpha_dbl / sigma_dbl;
75-
}
76-
}
77-
78-
if (!is_constant_all<T_y>::value) {
79-
for (size_t n = 0; n < stan::math::size(y); ++n) {
80-
ops_partials.edge1_.partials_[n] *= cdf;
81-
}
82-
}
83-
if (!is_constant_all<T_shape>::value) {
84-
for (size_t n = 0; n < stan::math::size(alpha); ++n) {
85-
ops_partials.edge2_.partials_[n] *= cdf;
86-
}
87-
}
88-
if (!is_constant_all<T_scale>::value) {
89-
for (size_t n = 0; n < stan::math::size(sigma); ++n) {
90-
ops_partials.edge3_.partials_[n] *= cdf;
94+
ops_partials.edge2_.partials_ = rep_deriv * log(y_val / sigma_val);
9195
}
9296
}
9397
return ops_partials.build(cdf);

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

+37-26
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@
66
#include <stan/math/prim/fun/log.hpp>
77
#include <stan/math/prim/fun/size_zero.hpp>
88
#include <stan/math/prim/fun/max_size.hpp>
9+
#include <stan/math/prim/fun/to_ref.hpp>
910
#include <stan/math/prim/fun/value_of.hpp>
1011
#include <stan/math/prim/functor/operands_and_partials.hpp>
1112
#include <cmath>
@@ -32,41 +33,51 @@ return_type_t<T_y, T_shape, T_scale> weibull_lccdf(const T_y& y,
3233
const T_shape& alpha,
3334
const T_scale& sigma) {
3435
using T_partials_return = partials_return_t<T_y, T_shape, T_scale>;
35-
using std::log;
36+
using T_y_ref = ref_type_if_t<!is_constant<T_y>::value, T_y>;
37+
using T_alpha_ref = ref_type_if_t<!is_constant<T_shape>::value, T_shape>;
38+
using T_sigma_ref = ref_type_if_t<!is_constant<T_scale>::value, T_scale>;
3639
using std::pow;
3740
static const char* function = "weibull_lccdf";
38-
check_nonnegative(function, "Random variable", y);
39-
check_positive_finite(function, "Shape parameter", alpha);
40-
check_positive_finite(function, "Scale parameter", sigma);
41+
42+
T_y_ref y_ref = y;
43+
T_alpha_ref alpha_ref = alpha;
44+
T_sigma_ref sigma_ref = sigma;
45+
46+
const auto& y_col = as_column_vector_or_scalar(y_ref);
47+
const auto& alpha_col = as_column_vector_or_scalar(alpha_ref);
48+
const auto& sigma_col = as_column_vector_or_scalar(sigma_ref);
49+
50+
const auto& y_arr = as_array_or_scalar(y_col);
51+
const auto& alpha_arr = as_array_or_scalar(alpha_col);
52+
const auto& sigma_arr = as_array_or_scalar(sigma_col);
53+
54+
ref_type_t<decltype(value_of(y_arr))> y_val = value_of(y_arr);
55+
ref_type_t<decltype(value_of(alpha_arr))> alpha_val = value_of(alpha_arr);
56+
ref_type_t<decltype(value_of(sigma_arr))> sigma_val = value_of(sigma_arr);
57+
58+
check_nonnegative(function, "Random variable", y_val);
59+
check_positive_finite(function, "Shape parameter", alpha_val);
60+
check_positive_finite(function, "Scale parameter", sigma_val);
4161

4262
if (size_zero(y, alpha, sigma)) {
4363
return 0.0;
4464
}
4565

46-
T_partials_return ccdf_log(0.0);
47-
operands_and_partials<T_y, T_shape, T_scale> ops_partials(y, alpha, sigma);
48-
49-
scalar_seq_view<T_y> y_vec(y);
50-
scalar_seq_view<T_scale> sigma_vec(sigma);
51-
scalar_seq_view<T_shape> alpha_vec(alpha);
52-
size_t N = max_size(y, sigma, alpha);
53-
for (size_t n = 0; n < N; n++) {
54-
const T_partials_return y_dbl = value_of(y_vec[n]);
55-
const T_partials_return sigma_dbl = value_of(sigma_vec[n]);
56-
const T_partials_return alpha_dbl = value_of(alpha_vec[n]);
57-
const T_partials_return pow_n = pow(y_dbl / sigma_dbl, alpha_dbl);
66+
operands_and_partials<T_y_ref, T_alpha_ref, T_sigma_ref> ops_partials(
67+
y_ref, alpha_ref, sigma_ref);
5868

59-
ccdf_log -= pow_n;
69+
const auto& pow_n = to_ref_if<!is_constant_all<T_y, T_shape, T_scale>::value>(
70+
pow(y_val / sigma_val, alpha_val));
71+
T_partials_return ccdf_log = -sum(pow_n);
6072

61-
if (!is_constant_all<T_y>::value) {
62-
ops_partials.edge1_.partials_[n] -= alpha_dbl / y_dbl * pow_n;
63-
}
64-
if (!is_constant_all<T_shape>::value) {
65-
ops_partials.edge2_.partials_[n] -= log(y_dbl / sigma_dbl) * pow_n;
66-
}
67-
if (!is_constant_all<T_scale>::value) {
68-
ops_partials.edge3_.partials_[n] += alpha_dbl / sigma_dbl * pow_n;
69-
}
73+
if (!is_constant_all<T_y>::value) {
74+
ops_partials.edge1_.partials_ = -alpha_val / y_val * pow_n;
75+
}
76+
if (!is_constant_all<T_shape>::value) {
77+
ops_partials.edge2_.partials_ = -log(y_val / sigma_val) * pow_n;
78+
}
79+
if (!is_constant_all<T_scale>::value) {
80+
ops_partials.edge3_.partials_ = alpha_val / sigma_val * pow_n;
7081
}
7182
return ops_partials.build(ccdf_log);
7283
}

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

+45-26
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77
#include <stan/math/prim/fun/log.hpp>
88
#include <stan/math/prim/fun/size_zero.hpp>
99
#include <stan/math/prim/fun/max_size.hpp>
10+
#include <stan/math/prim/fun/to_ref.hpp>
1011
#include <stan/math/prim/fun/value_of.hpp>
1112
#include <stan/math/prim/functor/operands_and_partials.hpp>
1213
#include <cmath>
@@ -33,44 +34,62 @@ return_type_t<T_y, T_shape, T_scale> weibull_lcdf(const T_y& y,
3334
const T_shape& alpha,
3435
const T_scale& sigma) {
3536
using T_partials_return = partials_return_t<T_y, T_shape, T_scale>;
36-
using std::exp;
37-
using std::log;
37+
using T_y_ref = ref_type_if_t<!is_constant<T_y>::value, T_y>;
38+
using T_alpha_ref = ref_type_if_t<!is_constant<T_shape>::value, T_shape>;
39+
using T_sigma_ref = ref_type_if_t<!is_constant<T_scale>::value, T_scale>;
3840
using std::pow;
3941
static const char* function = "weibull_lcdf";
40-
check_nonnegative(function, "Random variable", y);
41-
check_positive_finite(function, "Shape parameter", alpha);
42-
check_positive_finite(function, "Scale parameter", sigma);
42+
43+
T_y_ref y_ref = y;
44+
T_alpha_ref alpha_ref = alpha;
45+
T_sigma_ref sigma_ref = sigma;
46+
47+
const auto& y_col = as_column_vector_or_scalar(y_ref);
48+
const auto& alpha_col = as_column_vector_or_scalar(alpha_ref);
49+
const auto& sigma_col = as_column_vector_or_scalar(sigma_ref);
50+
51+
const auto& y_arr = as_array_or_scalar(y_col);
52+
const auto& alpha_arr = as_array_or_scalar(alpha_col);
53+
const auto& sigma_arr = as_array_or_scalar(sigma_col);
54+
55+
ref_type_t<decltype(value_of(y_arr))> y_val = value_of(y_arr);
56+
ref_type_t<decltype(value_of(alpha_arr))> alpha_val = value_of(alpha_arr);
57+
ref_type_t<decltype(value_of(sigma_arr))> sigma_val = value_of(sigma_arr);
58+
59+
check_nonnegative(function, "Random variable", y_val);
60+
check_positive_finite(function, "Shape parameter", alpha_val);
61+
check_positive_finite(function, "Scale parameter", sigma_val);
4362

4463
if (size_zero(y, alpha, sigma)) {
4564
return 0.0;
4665
}
4766

48-
T_partials_return cdf_log(0.0);
49-
operands_and_partials<T_y, T_shape, T_scale> ops_partials(y, alpha, sigma);
67+
operands_and_partials<T_y_ref, T_alpha_ref, T_sigma_ref> ops_partials(
68+
y_ref, alpha_ref, sigma_ref);
5069

51-
scalar_seq_view<T_y> y_vec(y);
52-
scalar_seq_view<T_scale> sigma_vec(sigma);
53-
scalar_seq_view<T_shape> alpha_vec(alpha);
54-
size_t N = max_size(y, sigma, alpha);
55-
for (size_t n = 0; n < N; n++) {
56-
const T_partials_return y_dbl = value_of(y_vec[n]);
57-
const T_partials_return sigma_dbl = value_of(sigma_vec[n]);
58-
const T_partials_return alpha_dbl = value_of(alpha_vec[n]);
59-
const T_partials_return pow_n = pow(y_dbl / sigma_dbl, alpha_dbl);
60-
const T_partials_return exp_n = exp(-pow_n);
61-
const T_partials_return cdf_n = 1.0 - exp_n;
70+
constexpr bool any_derivs = !is_constant_all<T_y, T_shape, T_scale>::value;
71+
const auto& pow_n = to_ref_if<any_derivs>(pow(y_val / sigma_val, alpha_val));
72+
const auto& exp_n = to_ref_if<any_derivs>(exp(-pow_n));
6273

63-
cdf_log += log(cdf_n);
74+
T_partials_return cdf_log = sum(log(1 - exp_n));
6475

65-
const T_partials_return rep_deriv = pow_n / (1.0 / exp_n - 1.0);
66-
if (!is_constant_all<T_y>::value) {
67-
ops_partials.edge1_.partials_[n] += rep_deriv * alpha_dbl / y_dbl;
76+
if (!is_constant_all<T_y, T_scale, T_shape>::value) {
77+
const auto& rep_deriv = to_ref_if<(!is_constant_all<T_y, T_scale>::value
78+
&& !is_constant_all<T_shape>::value)>(
79+
pow_n / (1.0 / exp_n - 1.0));
80+
if (!is_constant_all<T_y, T_scale>::value) {
81+
const auto& deriv_y_sigma = to_ref_if<(
82+
!is_constant_all<T_y>::value && !is_constant_all<T_scale>::value)>(
83+
rep_deriv * alpha_val);
84+
if (!is_constant_all<T_y>::value) {
85+
ops_partials.edge1_.partials_ = deriv_y_sigma / y_val;
86+
}
87+
if (!is_constant_all<T_scale>::value) {
88+
ops_partials.edge3_.partials_ = -deriv_y_sigma / sigma_val;
89+
}
6890
}
6991
if (!is_constant_all<T_shape>::value) {
70-
ops_partials.edge2_.partials_[n] += rep_deriv * log(y_dbl / sigma_dbl);
71-
}
72-
if (!is_constant_all<T_scale>::value) {
73-
ops_partials.edge3_.partials_[n] -= rep_deriv * alpha_dbl / sigma_dbl;
92+
ops_partials.edge2_.partials_ = rep_deriv * log(y_val / sigma_val);
7493
}
7594
}
7695
return ops_partials.build(cdf_log);

0 commit comments

Comments
 (0)