1 #ifndef STAN_MATH_PRIM_FUN_DIGAMMA_HPP
2 #define STAN_MATH_PRIM_FUN_DIGAMMA_HPP
3
4 #include <stan/math/prim/meta.hpp>
5 #include <stan/math/prim/fun/boost_policy.hpp>
6 #include <stan/math/prim/functor/apply_scalar_unary.hpp>
7 #include <boost/math/special_functions/digamma.hpp>
8
9 namespace stan {
10 namespace math {
11
12 /**
13 * Return the derivative of the log gamma function
14 * at the specified value.
15 *
16 \f[
17 \mbox{digamma}(x) =
18 \begin{cases}
19 \textrm{error} & \mbox{if } x\in \{\dots, -3, -2, -1, 0\}\\
20 \Psi(x) & \mbox{if } x\not\in \{\dots, -3, -2, -1, 0\}\\[6pt]
21 \textrm{NaN} & \mbox{if } x = \textrm{NaN}
22 \end{cases}
23 \f]
24
25 \f[
26 \frac{\partial\, \mbox{digamma}(x)}{\partial x} =
27 \begin{cases}
28 \textrm{error} & \mbox{if } x\in \{\dots, -3, -2, -1, 0\}\\
29 \frac{\partial\, \Psi(x)}{\partial x} & \mbox{if } x\not\in \{\dots, -3,
30 -2, -1, 0\}\\[6pt] \textrm{NaN} & \mbox{if } x = \textrm{NaN} \end{cases} \f]
31
32 \f[
33 \Psi(x)=\frac{\Gamma'(x)}{\Gamma(x)}
34 \f]
35
36 \f[
37 \frac{\partial \, \Psi(x)}{\partial x} =
38 \frac{\Gamma''(x)\Gamma(x)-(\Gamma'(x))^2}{\Gamma^2(x)} \f]
39
40 *
41 * The design follows the standard C++ library in returning NaN
42 * rather than throwing exceptions.
43 *
44 * @param[in] x argument
45 * @return derivative of log gamma function at argument
46 */
digamma(double x)47 inline double digamma(double x) {
48 return boost::math::digamma(x, boost_policy_t<>());
49 }
50
51 /**
52 * Structure to wrap digamma() so it can be vectorized.
53 *
54 * @tparam T type of variable
55 * @param x variable
56 * @return Digamma function applied to x.
57 * @throw std::domain_error if x is a negative integer or 0
58 */
59 struct digamma_fun {
60 template <typename T>
funstan::math::digamma_fun61 static inline T fun(const T& x) {
62 return digamma(x);
63 }
64 };
65
66 /**
67 * Vectorized version of digamma().
68 *
69 * @tparam T type of container
70 * @param x container
71 * @return Digamma function applied to each value in x.
72 * @throw std::domain_error if any value is a negative integer or 0
73 */
74 template <typename T,
75 require_not_nonscalar_prim_or_rev_kernel_expression_t<T>* = nullptr,
76 require_not_var_matrix_t<T>* = nullptr>
digamma(const T & x)77 inline auto digamma(const T& x) {
78 return apply_scalar_unary<digamma_fun, T>::apply(x);
79 }
80
81 } // namespace math
82 } // namespace stan
83
84 #endif
85