Stan Math Library  2.12.0
reverse mode automatic differentiation
trigamma.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_SCAL_FUN_TRIGAMMA_HPP
2 #define STAN_MATH_PRIM_SCAL_FUN_TRIGAMMA_HPP
3 
4  // Reference:
5  // BE Schneider,
6  // Algorithm AS 121:
7  // Trigamma Function,
8  // Applied Statistics,
9  // Volume 27, Number 1, pages 97-99, 1978.
10 
12 #include <cmath>
13 
14 namespace stan {
15  namespace math {
16 
46  template <typename T>
47  inline
48  T
49  trigamma(T x) {
50  using std::floor;
51  using std::sin;
52 
53  double small = 0.0001;
54  double large = 5.0;
55  T value;
56  T y;
57  T z;
58 
59  // bernoulli numbers
60  double b2 = 1.0 / 6.0;
61  double b4 = -1.0 / 30.0;
62  double b6 = 1.0 / 42.0;
63  double b8 = -1.0 / 30.0;
64 
65  // negative integers and zero return postiive infinity
66  // see http://mathworld.wolfram.com/PolygammaFunction.html
67  if ((x <= 0.0) && (floor(x) == x)) {
68  value = positive_infinity();
69  return value;
70  }
71 
72  // negative non-integers: use the reflection formula
73  // see http://mathworld.wolfram.com/PolygammaFunction.html
74  if ((x <= 0) && (floor(x) != x)) {
75  value = -trigamma(-x + 1.0) + (pi() / sin(-pi() * x))
76  * (pi() / sin(-pi() * x));
77  return value;
78  }
79 
80  // small value approximation if x <= small.
81  if (x <= small)
82  return 1.0 / (x * x);
83 
84  // use recurrence relation until x >= large
85  // see http://mathworld.wolfram.com/PolygammaFunction.html
86  z = x;
87  value = 0.0;
88  while (z < large) {
89  value += 1.0 / (z * z);
90  z += 1.0;
91  }
92 
93  // asymptotic expansion as a Laurent series if x >= large
94  // see http://en.wikipedia.org/wiki/Trigamma_function
95  y = 1.0 / (z * z);
96  value += 0.5 * y + (1.0 + y * (b2 + y * (b4 + y * (b6 + y * b8)))) / z;
97 
98  return value;
99  }
100  }
101 }
102 
103 #endif
T trigamma(T x)
Definition: trigamma.hpp:49
fvar< T > sin(const fvar< T > &x)
Definition: sin.hpp:12
double positive_infinity()
Return positive infinity.
Definition: constants.hpp:121
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:11
double pi()
Return the value of pi.
Definition: constants.hpp:85

     [ Stan Home Page ] © 2011–2016, Stan Development Team.