Stan Math Library  2.9.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 
16  namespace math {
17 
47  template <typename T>
48  inline
49  T
50  trigamma(T x) {
51  using std::floor;
52  using std::sin;
53 
54  double small = 0.0001;
55  double large = 5.0;
56  T value;
57  T y;
58  T z;
59 
60  // bernoulli numbers
61  double b2 = 1.0 / 6.0;
62  double b4 = -1.0 / 30.0;
63  double b6 = 1.0 / 42.0;
64  double b8 = -1.0 / 30.0;
65 
66  // negative integers and zero return postiive infinity
67  // see http:// mathworld.wolfram.com/PolygammaFunction.html
68  if ((x <= 0.0) && (floor(x) == x)) {
69  value = positive_infinity();
70  return value;
71  }
72 
73  // negative non-integers: use the reflection formula
74  // see http:// mathworld.wolfram.com/PolygammaFunction.html
75  if ((x <= 0) && (floor(x) != x)) {
76  value = -trigamma(-x + 1.0) + (pi() / sin(-pi() * x))
77  * (pi() / sin(-pi() * x));
78  return value;
79  }
80 
81  // small value approximation if x <= small.
82  if (x <= small)
83  return 1.0 / (x * x);
84 
85  // use recurrence relation until x >= large
86  // see http:// mathworld.wolfram.com/PolygammaFunction.html
87  z = x;
88  value = 0.0;
89  while (z < large) {
90  value += 1.0 / (z * z);
91  z += 1.0;
92  }
93 
94  // asymptotic expansion as a Laurent series if x >= large
95  // see http:// en.wikipedia.org/wiki/Trigamma_function
96  y = 1.0 / (z * z);
97  value += 0.5 * y + (1.0 + y * (b2 + y * (b4 + y * (b6 + y * b8)))) / z;
98 
99  return value;
100  }
101  }
102 }
103 
104 #endif
T trigamma(T x)
Definition: trigamma.hpp:50
fvar< T > sin(const fvar< T > &x)
Definition: sin.hpp:14
double positive_infinity()
Return positive infinity.
Definition: constants.hpp:123
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:11
double pi()
Return the value of pi.
Definition: constants.hpp:86

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