Stan Math Library  2.20.0
reverse mode automatic differentiation
log_modified_bessel_first_kind.hpp
Go to the documentation of this file.
1 // Copyright (c) 2006 Xiaogang Zhang
2 // Copyright (c) 2007, 2017 John Maddock
3 
5 #ifndef STAN_MATH_PRIM_SCAL_FUN_LOG_MODIFIED_BESSEL_FIRST_KIND_HPP
6 #define STAN_MATH_PRIM_SCAL_FUN_LOG_MODIFIED_BESSEL_FIRST_KIND_HPP
7 
8 #include <boost/math/tools/rational.hpp>
20 
21 namespace stan {
22 namespace math {
23 
24 /* Log of the modified Bessel function of the first kind,
25  * which is better known as the Bessel I function. See
26  * modified_bessel_first_kind.hpp for the function definition.
27  * The derivatives are known to be incorrect for v = 0 because a
28  * simple constant 0 is returned.
29  *
30  * @tparam T1 type of the order (v)
31  * @tparam T2 type of argument (z)
32  * @param v Order, can be a non-integer but must be at least -1
33  * @param z Real non-negative number
34  * @throws std::domain_error if either v or z is NaN, z is
35  * negative, or v is less than -1
36  * @return log of Bessel I function
37  */
38 template <typename T1, typename T2>
39 inline typename boost::math::tools::promote_args<T1, T2, double>::type
40 log_modified_bessel_first_kind(const T1 v, const T2 z) {
41  check_not_nan("log_modified_bessel_first_kind", "first argument (order)", v);
42  check_not_nan("log_modified_bessel_first_kind", "second argument (variable)",
43  z);
44  check_nonnegative("log_modified_bessel_first_kind",
45  "second argument (variable)", z);
46  check_greater_or_equal("log_modified_bessel_first_kind",
47  "first argument (order)", v, -1);
48 
49  using boost::math::tools::evaluate_polynomial;
50  using std::log;
51  using std::sqrt;
52 
53  typedef typename boost::math::tools::promote_args<T1, T2, double>::type T;
54 
55  if (z == 0) {
56  if (v == 0)
57  return 0.0;
58  if (v > 0)
59  return NEGATIVE_INFTY;
60  return INFTY;
61  }
62  if (is_inf(z))
63  return z;
64  if (v == 0) {
65  // WARNING: will not autodiff for v = 0 correctly
66  // modified from Boost's bessel_i0_imp in the double precision case,
67  // which refers to:
68  // Modified Bessel function of the first kind of order zero
69  // we use the approximating forms derived in:
70  // "Rational Approximations for the Modified Bessel Function of the
71  // First Kind -- I0(x) for Computations with Double Precision"
72  // by Pavel Holoborodko, see
73  // http://www.advanpix.com/2015/11/11/rational-approximations-for-the-modified-bessel-function-of-the-first-kind-i0-computations-double-precision
74  // The actual coefficients used are [Boost's] own, and extend
75  // Pavel's work to precisions other than double.
76 
77  if (z < 7.75) {
78  // Bessel I0 over[10 ^ -16, 7.75]
79  // Max error in interpolated form : 3.042e-18
80  // Max Error found at double precision = Poly : 5.106609e-16
81  // Cheb : 5.239199e-16
82  static const double P[]
83  = {1.00000000000000000e+00, 2.49999999999999909e-01,
84  2.77777777777782257e-02, 1.73611111111023792e-03,
85  6.94444444453352521e-05, 1.92901234513219920e-06,
86  3.93675991102510739e-08, 6.15118672704439289e-10,
87  7.59407002058973446e-12, 7.59389793369836367e-14,
88  6.27767773636292611e-16, 4.34709704153272287e-18,
89  2.63417742690109154e-20, 1.13943037744822825e-22,
90  9.07926920085624812e-25};
91  return log1p_exp(multiply_log(2, z) - log(4.0)
92  + log(evaluate_polynomial(P, 0.25 * square(z))));
93  }
94  if (z < 500) {
95  // Max error in interpolated form : 1.685e-16
96  // Max Error found at double precision = Poly : 2.575063e-16
97  // Cheb : 2.247615e+00
98  static const double P[]
99  = {3.98942280401425088e-01, 4.98677850604961985e-02,
100  2.80506233928312623e-02, 2.92211225166047873e-02,
101  4.44207299493659561e-02, 1.30970574605856719e-01,
102  -3.35052280231727022e+00, 2.33025711583514727e+02,
103  -1.13366350697172355e+04, 4.24057674317867331e+05,
104  -1.23157028595698731e+07, 2.80231938155267516e+08,
105  -5.01883999713777929e+09, 7.08029243015109113e+10,
106  -7.84261082124811106e+11, 6.76825737854096565e+12,
107  -4.49034849696138065e+13, 2.24155239966958995e+14,
108  -8.13426467865659318e+14, 2.02391097391687777e+15,
109  -3.08675715295370878e+15, 2.17587543863819074e+15};
110  return z + log(evaluate_polynomial(P, inv(z))) - multiply_log(0.5, z);
111  }
112  // Max error in interpolated form : 2.437e-18
113  // Max Error found at double precision = Poly : 1.216719e-16
114  static const double P[] = {3.98942280401432905e-01, 4.98677850491434560e-02,
115  2.80506308916506102e-02, 2.92179096853915176e-02,
116  4.53371208762579442e-02};
117  return z + log(evaluate_polynomial(P, inv(z))) - multiply_log(0.5, z);
118  }
119  if (v == 1) { // WARNING: will not autodiff for v = 1 correctly
120  // modified from Boost's bessel_i1_imp in the double precision case
121  // see credits above in the v == 0 case
122  if (z < 7.75) {
123  // Bessel I0 over[10 ^ -16, 7.75]
124  // Max error in interpolated form: 5.639e-17
125  // Max Error found at double precision = Poly: 1.795559e-16
126 
127  static const double P[]
128  = {8.333333333333333803e-02, 6.944444444444341983e-03,
129  3.472222222225921045e-04, 1.157407407354987232e-05,
130  2.755731926254790268e-07, 4.920949692800671435e-09,
131  6.834657311305621830e-11, 7.593969849687574339e-13,
132  6.904822652741917551e-15, 5.220157095351373194e-17,
133  3.410720494727771276e-19, 1.625212890947171108e-21,
134  1.332898928162290861e-23};
135  T a = square(z) * 0.25;
136  T Q[3] = {1, 0.5, evaluate_polynomial(P, a)};
137  return log(z) + log(evaluate_polynomial(Q, a)) - LOG_2;
138  }
139  if (z < 500) {
140  // Max error in interpolated form: 1.796e-16
141  // Max Error found at double precision = Poly: 2.898731e-16
142 
143  static const double P[]
144  = {3.989422804014406054e-01, -1.496033551613111533e-01,
145  -4.675104253598537322e-02, -4.090895951581637791e-02,
146  -5.719036414430205390e-02, -1.528189554374492735e-01,
147  3.458284470977172076e+00, -2.426181371595021021e+02,
148  1.178785865993440669e+04, -4.404655582443487334e+05,
149  1.277677779341446497e+07, -2.903390398236656519e+08,
150  5.192386898222206474e+09, -7.313784438967834057e+10,
151  8.087824484994859552e+11, -6.967602516005787001e+12,
152  4.614040809616582764e+13, -2.298849639457172489e+14,
153  8.325554073334618015e+14, -2.067285045778906105e+15,
154  3.146401654361325073e+15, -2.213318202179221945e+15};
155  return z + log(evaluate_polynomial(P, inv(z))) - multiply_log(0.5, z);
156  }
157  // Max error in interpolated form: 1.320e-19
158  // Max Error found at double precision = Poly: 7.065357e-17
159  static const double P[]
160  = {3.989422804014314820e-01, -1.496033551467584157e-01,
161  -4.675105322571775911e-02, -4.090421597376992892e-02,
162  -5.843630344778927582e-02};
163  return z + log(evaluate_polynomial(P, inv(z))) - multiply_log(0.5, z);
164  }
165  if (z > 100) {
166  // Boost does something like this in asymptotic_bessel_i_large_x
167  T lim = pow((square(v) + 2.5) / (2 * z), 3) / 24;
168  if (lim < (EPSILON * 10)) {
169  T s = 1;
170  T mu = 4 * square(v);
171  T ex = 8 * z;
172  T num = mu - 1;
173  T denom = ex;
174  s -= num / denom;
175  num *= mu - 9;
176  denom *= ex * 2;
177  s += num / denom;
178  num *= mu - 25;
179  denom *= ex * 3;
180  s -= num / denom;
181  s = z - log(sqrt(2 * z * pi())) + log(s);
182  return s;
183  }
184  }
185 
186  typename boost::math::tools::promote_args<T2>::type log_half_z = log(0.5 * z);
187  typename boost::math::tools::promote_args<T1>::type lgam
188  = v > -1 ? lgamma(v + 1.0) : 0;
189  T lcons = (2.0 + v) * log_half_z;
190  T out;
191  if (v > -1) {
192  out = log_sum_exp(v * log_half_z - lgam, lcons - lgamma(v + 2));
193  lgam += log1p(v);
194  } else {
195  out = lcons;
196  }
197 
198  double m = 2;
199  double lfac = 0;
200  T old_out;
201  do {
202  lfac += log(m);
203  lgam += log(v + m);
204  lcons += 2 * log_half_z;
205  old_out = out;
206  out = log_sum_exp(out, lcons - lfac - lgam); // underflows eventually
207  m++;
208  } while (out > old_out || out < old_out);
209  return out;
210 }
211 
212 } // namespace math
213 } // namespace stan
214 
215 #endif
const double LOG_2
The natural logarithm of 2, .
Definition: constants.hpp:37
fvar< T > lgamma(const fvar< T > &x)
Return the natural logarithm of the gamma function applied to the specified argument.
Definition: lgamma.hpp:21
fvar< T > sqrt(const fvar< T > &x)
Definition: sqrt.hpp:13
fvar< T > log(const fvar< T > &x)
Definition: log.hpp:12
fvar< T > log_sum_exp(const std::vector< fvar< T > > &v)
Definition: log_sum_exp.hpp:12
void check_nonnegative(const char *function, const char *name, const T_y &y)
Check if y is non-negative.
void check_greater_or_equal(const char *function, const char *name, const T_y &y, const T_low &low)
Check if y is greater or equal than low.
fvar< T > square(const fvar< T > &x)
Definition: square.hpp:12
const double EPSILON
Smallest positive value.
Definition: constants.hpp:63
void check_not_nan(const char *function, const char *name, const T_y &y)
Check if y is not NaN.
int is_inf(const fvar< T > &x)
Returns 1 if the input&#39;s value is infinite and 0 otherwise.
Definition: is_inf.hpp:20
fvar< T > multiply_log(const fvar< T > &x1, const fvar< T > &x2)
fvar< T > log1p_exp(const fvar< T > &x)
Definition: log1p_exp.hpp:12
boost::math::tools::promote_args< T1, T2, double >::type log_modified_bessel_first_kind(const T1 v, const T2 z)
fvar< T > log1p(const fvar< T > &x)
Definition: log1p.hpp:12
const double INFTY
Positive infinity.
Definition: constants.hpp:48
double pi()
Return the value of pi.
Definition: constants.hpp:80
fvar< T > pow(const fvar< T > &x1, const fvar< T > &x2)
Definition: pow.hpp:16
const double NEGATIVE_INFTY
Negative infinity.
Definition: constants.hpp:53
fvar< T > inv(const fvar< T > &x)
Definition: inv.hpp:12

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