| Stan Math Library 2.12.0 reverse mode automatic differentiation

- Main Page
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

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
11 #include <stan/math/prim/scal/fun/constants.hpp>
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
stan::math::trigamma
T trigamma(T x)
Definition: trigamma.hpp:49
stan
Definition: log_sum_exp.hpp:8
stan::math::sin
fvar< T > sin(const fvar< T > &x)
Definition: sin.hpp:12
stan::math::positive_infinity
double positive_infinity()
Return positive infinity.
Definition: constants.hpp:121
constants.hpp
stan::math::floor
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:11
stan::math::pi
double pi()
Return the value of pi.
Definition: constants.hpp:85

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