1 #ifndef STAN_MATH_PRIM_FUN_LMGAMMA_HPP
2 #define STAN_MATH_PRIM_FUN_LMGAMMA_HPP
3
4 #include <stan/math/prim/meta.hpp>
5 #include <stan/math/prim/fun/constants.hpp>
6 #include <stan/math/prim/fun/lgamma.hpp>
7 #include <stan/math/prim/fun/sum.hpp>
8 #include <stan/math/prim/functor/apply_scalar_binary.hpp>
9
10 namespace stan {
11 namespace math {
12
13 /**
14 * Return the natural logarithm of the multivariate gamma function
15 * with the specified dimensions and argument.
16 *
17 * <p>The multivariate gamma function \f$\Gamma_k(x)\f$ for
18 * dimensionality \f$k\f$ and argument \f$x\f$ is defined by
19 *
20 * <p>\f$\Gamma_k(x) = \pi^{k(k-1)/4} \, \prod_{j=1}^k \Gamma(x + (1 - j)/2)\f$,
21 *
22 * where \f$\Gamma()\f$ is the gamma function.
23 *
24 \f[
25 \mbox{lmgamma}(n, x) =
26 \begin{cases}
27 \textrm{error} & \mbox{if } x\in \{\dots, -3, -2, -1, 0\}\\
28 \ln\Gamma_n(x) & \mbox{if } x\not\in \{\dots, -3, -2, -1, 0\}\\[6pt]
29 \textrm{NaN} & \mbox{if } x = \textrm{NaN}
30 \end{cases}
31 \f]
32
33 \f[
34 \frac{\partial\, \mbox{lmgamma}(n, x)}{\partial x} =
35 \begin{cases}
36 \textrm{error} & \mbox{if } x\in \{\dots, -3, -2, -1, 0\}\\
37 \frac{\partial\, \ln\Gamma_n(x)}{\partial x} & \mbox{if } x\not\in \{\dots,
38 -3, -2, -1, 0\}\\[6pt] \textrm{NaN} & \mbox{if } x = \textrm{NaN} \end{cases}
39 \f]
40
41 \f[
42 \ln\Gamma_n(x) = \pi^{n(n-1)/4} \, \prod_{j=1}^n \Gamma(x + (1 - j)/2)
43 \f]
44
45 \f[
46 \frac{\partial \, \ln\Gamma_n(x)}{\partial x} = \sum_{j=1}^n \Psi(x + (1 - j)
47 / 2) \f]
48 *
49 * @tparam T type of scalar
50 * @param k Number of dimensions.
51 * @param x Function argument.
52 * @return Natural log of the multivariate gamma function.
53 */
54 template <typename T, require_arithmetic_t<T>* = nullptr>
lmgamma(int k,T x)55 inline return_type_t<T> lmgamma(int k, T x) {
56 return_type_t<T> result = k * (k - 1) * LOG_PI_OVER_FOUR;
57
58 return result + sum(lgamma(x + (1 - Eigen::ArrayXd::LinSpaced(k, 1, k)) / 2));
59 }
60
61 /**
62 * Enables the vectorised application of the natural log of the multivariate
63 * gamma function, when the first and/or second arguments are containers.
64 *
65 * @tparam T1 type of first input
66 * @tparam T2 type of second input
67 * @param a First input
68 * @param b Second input
69 * @return Natural log of the multivariate gamma function applied to the two
70 * inputs.
71 */
72 template <typename T1, typename T2, require_any_container_t<T1, T2>* = nullptr>
lmgamma(const T1 & a,const T2 & b)73 inline auto lmgamma(const T1& a, const T2& b) {
74 return apply_scalar_binary(
75 a, b, [&](const auto& c, const auto& d) { return lmgamma(c, d); });
76 }
77
78 } // namespace math
79 } // namespace stan
80 #endif
81