1 #ifndef STAN_MATH_PRIM_FUN_LGAMMA_STIRLING_DIFF_HPP
2 #define STAN_MATH_PRIM_FUN_LGAMMA_STIRLING_DIFF_HPP
3
4 #include <stan/math/prim/meta.hpp>
5 #include <stan/math/prim/err.hpp>
6 #include <stan/math/prim/fun/constants.hpp>
7 #include <stan/math/prim/fun/inv.hpp>
8 #include <stan/math/prim/fun/is_nan.hpp>
9 #include <stan/math/prim/fun/lgamma.hpp>
10 #include <stan/math/prim/fun/lgamma_stirling.hpp>
11 #include <stan/math/prim/fun/square.hpp>
12 #include <stan/math/prim/fun/value_of.hpp>
13 #include <cmath>
14
15 namespace stan {
16 namespace math {
17
18 constexpr double lgamma_stirling_diff_useful = 10;
19
20 /**
21 * Return the difference between log of the gamma function and its Stirling
22 * approximation.
23 * This is useful to stably compute log of ratios of gamma functions with large
24 * arguments where the Stirling approximation allows for analytic solution
25 * and the (small) differences can be added afterwards.
26 * This is for example used in the implementation of lbeta.
27 *
28 * The function will return correct value for all arguments, but using it can
29 * lead to a loss of precision when x < lgamma_stirling_diff_useful.
30 *
31 \f[
32 \mbox{lgamma_stirling_diff}(x) =
33 \log(\Gamma(x)) - \frac{1}{2} \log(2\pi) +
34 (x-\frac{1}{2})*\log(x) - x
35 \f]
36
37 *
38 * @tparam T type of value
39 * @param x value
40 * @return Difference between lgamma(x) and its Stirling approximation.
41 */
42 template <typename T>
lgamma_stirling_diff(const T x)43 return_type_t<T> lgamma_stirling_diff(const T x) {
44 using T_ret = return_type_t<T>;
45
46 if (is_nan(value_of_rec(x))) {
47 return NOT_A_NUMBER;
48 }
49 check_nonnegative("lgamma_stirling_diff", "argument", x);
50
51 if (x == 0) {
52 return INFTY;
53 }
54 if (value_of(x) < lgamma_stirling_diff_useful) {
55 return lgamma(x) - lgamma_stirling(x);
56 }
57
58 // Using the Stirling series as expressed in formula 5.11.1. at
59 // https://dlmf.nist.gov/5.11
60 constexpr double stirling_series[]{
61 0.0833333333333333333333333, -0.00277777777777777777777778,
62 0.000793650793650793650793651, -0.000595238095238095238095238,
63 0.000841750841750841750841751, -0.00191752691752691752691753,
64 0.00641025641025641025641026, -0.0295506535947712418300654};
65
66 constexpr int n_stirling_terms = 6;
67 T_ret result(0.0);
68 T_ret multiplier = inv(x);
69 T_ret inv_x_squared = square(multiplier);
70 for (int n = 0; n < n_stirling_terms; n++) {
71 if (n > 0) {
72 multiplier *= inv_x_squared;
73 }
74 result += stirling_series[n] * multiplier;
75 }
76 return result;
77 }
78
79 } // namespace math
80 } // namespace stan
81
82 #endif
83