TMB Documentation  v1.9.11
tweedie.cpp
1 // Re-structured version of tweedie density function from 'cplm' package.
2
3 /* the threshold used in finding the bounds of the series */
4 #define TWEEDIE_DROP 37.0
5 /* the loop increment used in finding the bounds of the series */
6 #define TWEEDIE_INCRE 5
7 /* the max number of terms allowed in the finite sum approximation*/
8 #define TWEEDIE_NTERM 20000
9
18 template<class Float>
19 Float tweedie_logW(Float y, Float phi, Float p){
20  bool ok = (0 < y) && (0 < phi) && (1 < p) && (p < 2);
21  if (!ok) return NAN;
22
23  Float p1 = p - 1.0, p2 = 2.0 - p;
24  Float a = - p2 / p1, a1 = 1.0 / p1;
25  Float cc, w, sum_ww = 0.0;
26  double ww_max = -INFINITY ;
27  double j;
28
29  /* only need the lower bound and the # terms to be stored */
30  double jh, jl, jd;
31  double jmax = 0;
32  Float logz = 0;
33
34  /* compute jmax for the given y > 0*/
35  cc = a * log(p1) - log(p2);
36  jmax = asDouble( fmax2(1.0, pow(y, p2) / (phi * p2)) );
37  logz = - a * log(y) - a1 * log(phi) + cc;
38
39  /* find bounds in the summation */
40  /* locate upper bound */
41  cc = logz + a1 + a * log(-a);
42  j = jmax ;
43  w = a1 * j ;
44  while (1) {
45  j += TWEEDIE_INCRE ;
46  if (j * (cc - a1 * log(j)) < (w - TWEEDIE_DROP))
47  break ;
48  }
49  jh = ceil(j);
50  /* locate lower bound */
51  j = jmax;
52  while (1) {
53  j -= TWEEDIE_INCRE ;
54  if (j < 1 || j * (cc - a1 * log(j)) < w - TWEEDIE_DROP)
55  break ;
56  }
57  jl = std::fmax(1., floor(j)) ;
58  jd = jh - jl + 1;
59
60  /* set limit for # terms in the sum */
61  int nterms = (int) std::fmin(jd, TWEEDIE_NTERM) ;
62  std::vector<Float> ww(nterms);
63  /* evaluate series using the finite sum*/
64  /* y > 0 */
65  sum_ww = 0.0 ;
66  int iterm = (int) std::fmin(jd, nterms) ;
67  for (int k = 0; k < iterm; k++) {
68  j = k + jl ;
69  ww[k] = j * logz - lgamma(1 + j) - lgamma(-a * j);
70  ww_max = std::fmax( ww_max, asDouble(ww[k]) );
71  }
72  for (int k = 0; k < iterm; k++)
73  sum_ww += exp(ww[k] - ww_max);
74  Float ans = log(sum_ww) + ww_max ;
75
76  return ans;
77 }
Type lgamma(Type x)
Logarithm of gamma function (following R argument convention).
Definition: lgamma.hpp:12