Automatic Differentiation
 
Loading...
Searching...
No Matches
lbeta.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_FUN_LBETA_HPP
2#define STAN_MATH_PRIM_FUN_LBETA_HPP
3
14
15namespace stan {
16namespace math {
17
64template <typename T1, typename T2, require_all_arithmetic_t<T1, T2>* = nullptr>
65return_type_t<T1, T2> lbeta(const T1 a, const T2 b) {
66 using T_ret = return_type_t<T1, T2>;
67
68 if (is_any_nan(a, b)) {
69 return NOT_A_NUMBER;
70 }
71
72 static constexpr const char* function = "lbeta";
73 check_nonnegative(function, "first argument", a);
74 check_nonnegative(function, "second argument", b);
75 T_ret x; // x is the smaller of the two
76 T_ret y;
77 if (a < b) {
78 x = a;
79 y = b;
80 } else {
81 x = b;
82 y = a;
83 }
84
85 // Special cases
86 if (x == 0) {
87 return INFTY;
88 }
89 if (is_inf(y)) {
90 return NEGATIVE_INFTY;
91 }
92
93 // For large x or y, separate the lgamma values into Stirling approximations
94 // and appropriate corrections. The Stirling approximations allow for
95 // analytic simplification and the corrections are added later.
96 //
97 // The overall approach is inspired by the code in R, where the algorithm is
98 // credited to W. Fullerton of Los Alamos Scientific Laboratory
100 // both small
101 return lgamma(x) + lgamma(y) - lgamma(x + y);
102 }
103 T_ret x_over_xy = x / (x + y);
105 // y large, x small
106 T_ret stirling_diff = lgamma_stirling_diff(y) - lgamma_stirling_diff(x + y);
107 T_ret stirling = (y - 0.5) * log1m(x_over_xy) + x * (1 - log(x + y));
108 return stirling + lgamma(x) + stirling_diff;
109 }
110
111 // both large
112 T_ret stirling_diff = lgamma_stirling_diff(x) + lgamma_stirling_diff(y)
113 - lgamma_stirling_diff(x + y);
114 T_ret stirling = (x - 0.5) * log(x_over_xy) + y * log1m(x_over_xy)
115 + HALF_LOG_TWO_PI - 0.5 * log(y);
116 return stirling + stirling_diff;
117}
118
129template <typename T1, typename T2, require_any_container_t<T1, T2>* = nullptr>
130inline auto lbeta(const T1& a, const T2& b) {
131 return apply_scalar_binary(
132 a, b, [&](const auto& c, const auto& d) { return lbeta(c, d); });
133}
134
135} // namespace math
136} // namespace stan
137
138#endif
typename return_type< Ts... >::type return_type_t
Convenience type for the return type of the specified template parameters.
static constexpr double NOT_A_NUMBER
(Quiet) not-a-number value.
Definition constants.hpp:56
void check_nonnegative(const char *function, const char *name, const T_y &y)
Check if y is non-negative.
fvar< T > lbeta(const fvar< T > &x1, const fvar< T > &x2)
Definition lbeta.hpp:14
fvar< T > log(const fvar< T > &x)
Definition log.hpp:15
static constexpr double NEGATIVE_INFTY
Negative infinity.
Definition constants.hpp:51
fvar< T > lgamma(const fvar< T > &x)
Return the natural logarithm of the gamma function applied to the specified argument.
Definition lgamma.hpp:21
return_type_t< T > lgamma_stirling_diff(const T x)
Return the difference between log of the gamma function and its Stirling approximation.
int is_inf(const fvar< T > &x)
Returns 1 if the input's value is infinite and 0 otherwise.
Definition is_inf.hpp:21
constexpr double lgamma_stirling_diff_useful
static constexpr double HALF_LOG_TWO_PI
The value of half the natural logarithm , .
auto apply_scalar_binary(const T1 &x, const T2 &y, const F &f)
Base template function for vectorization of binary scalar functions defined by applying a functor to ...
fvar< T > log1m(const fvar< T > &x)
Definition log1m.hpp:12
bool is_any_nan(const T &x)
Returns true if the input is NaN and false otherwise.
static constexpr double INFTY
Positive infinity.
Definition constants.hpp:46
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...