TMB Documentation  v1.9.1
lgamma.hpp
Go to the documentation of this file.
1 // Copyright (C) 2013-2015 Kasper Kristensen
2 // License: GPL-2
3 
11 template<class Type>
12 Type lgamma(Type x){
13  CppAD::vector<Type> tx(2);
14  tx[0] = x;
15  tx[1] = Type(0);
16  return atomic::D_lgamma(tx)[0];
17 }
19 
20 
23 template<class Type>
24 Type lfactorial(Type x){
25  CppAD::vector<Type> tx(2);
26  tx[0] = x + Type(1);
27  tx[1] = Type(0);
28  return atomic::D_lgamma(tx)[0];
29 }
31 
32 /* Old lgamma approximation */
33 template <class Type>
34 inline Type lgamma_approx(const Type &y)
35 {
36  /* coefficients for gamma=7, kmax=8 Lanczos method */
37  static const Type
38  LogRootTwoPi_ = 0.9189385332046727418,
39  lanczos_7_c[9] = {
40  0.99999999999980993227684700473478,
41  676.520368121885098567009190444019,
42  -1259.13921672240287047156078755283,
43  771.3234287776530788486528258894,
44  -176.61502916214059906584551354,
45  12.507343278686904814458936853,
46  -0.13857109526572011689554707,
47  9.984369578019570859563e-6,
48  1.50563273514931155834e-7
49  };
50  Type x=y;
51  int k;
52  Type Ag;
53  Type term1, term2;
54  x -= Type(1.0); /* Lanczos writes z! instead of Gamma(z) */
55  Ag = lanczos_7_c[0];
56  for(k=1; k<=8; k++) { Ag += lanczos_7_c[k]/(x+k); }
57  /* (x+0.5)*log(x+7.5) - (x+7.5) + LogRootTwoPi_ + log(Ag(x)) */
58  term1 = (x+Type(0.5))*log((x+Type(7.5))/Type(M_E));
59  term2 = LogRootTwoPi_ + log(Ag);
60  return term1 + (term2 - Type(7.0));
61 }
62 
68 template<class Type>
69 inline Type dnbinom(const Type &x, const Type &size, const Type &prob,
70  int give_log=0)
71 {
72  Type n=size;
73  Type p=prob;
74  Type logres = lgamma(x+n)-lgamma(n)-lgamma(x+Type(1))+
75  n*log(p)+x*log(Type(1)-p);
76  if (give_log) return logres; else return exp(logres);
77 }
79 
80 
85 template<class Type>
86 inline Type dnbinom2(const Type &x, const Type &mu, const Type &var,
87  int give_log=0)
88 {
89  Type p=mu/var;
90  Type n=mu*p/(Type(1)-p);
91  return dnbinom(x,n,p,give_log);
92 }
94 
95 
102 template<class Type>
103 inline Type dnbinom_robust(const Type &x,
104  const Type &log_mu,
105  const Type &log_var_minus_mu,
106  int give_log=0)
107 {
108  CppAD::vector<Type> tx(4);
109  tx[0] = x;
110  tx[1] = log_mu;
111  tx[2] = log_var_minus_mu;
112  tx[3] = 0;
113  Type ans = atomic::log_dnbinom_robust(tx)[0];
114  return ( give_log ? ans : exp(ans) );
115 }
117 
118 
122 template<class Type>
123 inline Type dpois(const Type &x, const Type &lambda, int give_log=0)
124 {
125  Type logres = -lambda + x*log(lambda) - lgamma(x+Type(1));
126  if (give_log) return logres; else return exp(logres);
127 }
129 
130 
133 template<class Type>
134 Type dgamma(Type y, Type shape, Type scale, int give_log=0)
135 {
136  Type logres=-lgamma(shape)+(shape-Type(1.0))*log(y)-y/scale-shape*log(scale);
137  if(give_log)return logres; else return exp(logres);
138 }
140 
141 
144 template<class Type>
145 inline Type dlgamma(Type y, Type shape, Type scale, int give_log=0)
146 {
147  Type logres=-lgamma(shape)-shape*log(scale)-exp(y)/scale+shape*y;
148  if(give_log)return logres; else return exp(logres);
149 }
151 
152 
157 template<class Type>
158 inline Type dzipois(const Type &x, const Type &lambda, const Type &zip, int give_log=0)
159 {
160  Type logres;
161  if (x==Type(0)) logres=log(zip + (Type(1)-zip)*dpois(x, lambda, false));
162  else logres=log(Type(1)-zip) + dpois(x, lambda, true);
163  if (give_log) return logres; else return exp(logres);
164 }
166 
167 
174 template<class Type>
175 inline Type dzinbinom(const Type &x, const Type &size, const Type &p, const Type & zip,
176  int give_log=0)
177 {
178  Type logres;
179  if (x==Type(0)) logres=log(zip + (Type(1)-zip)*dnbinom(x, size, p, false));
180  else logres=log(Type(1)-zip) + dnbinom(x, size, p, true);
181  if (give_log) return logres; else return exp(logres);
182 }
183 
191 template<class Type>
192 inline Type dzinbinom2(const Type &x, const Type &mu, const Type &var, const Type & zip,
193  int give_log=0)
194 {
195  Type p=mu/var;
196  Type n=mu*p/(Type(1)-p);
197  return dzinbinom(x,n,p,zip,give_log);
198 }
199 
200 /********************************************************************/
201 /* SIMULATON CODE */
202 /********************************************************************/
203 
204 extern "C" {
205  double Rf_rnbinom(double n, double p);
206 }
208 template<class Type>
209 Type rnbinom(Type n, Type p)
210 {
211  return Rf_rnbinom(asDouble(n), asDouble(p));
212 }
215 
216 
217 template<class Type>
218 Type rnbinom2(Type mu, Type var)
219 {
220  Type p = mu / var;
221  Type n = mu * p / (Type(1) - p);
222  return Rf_rnbinom(asDouble(n), asDouble(p));
223 }
#define VECTORIZE2_tt(FUN)
Vectorize 2-argument functions.
Definition: Vectorize.hpp:71
Type rnbinom(Type n, Type p)
Simulate from a negative binomial distribution.
Definition: lgamma.hpp:209
#define VECTORIZE2_n(FUN)
Add the &#39;n&#39; integer argument to a simulation method with two arguments.
Definition: Vectorize.hpp:154
Type dpois(const Type &x, const Type &lambda, int give_log=0)
Poisson probability function.
Definition: lgamma.hpp:123
#define VECTORIZE1_t(FUN)
Vectorize 1-argument functions.
Definition: Vectorize.hpp:63
Type dlgamma(Type y, Type shape, Type scale, int give_log=0)
Density of log(X) where X~gamma distributed.
Definition: lgamma.hpp:145
Type lfactorial(Type x)
Logarithm of factorial function (following R argument convention).
Definition: lgamma.hpp:24
Type dgamma(Type y, Type shape, Type scale, int give_log=0)
Density of X where X~gamma distributed.
Definition: lgamma.hpp:134
Type dzinbinom2(const Type &x, const Type &mu, const Type &var, const Type &zip, int give_log=0)
Zero-Inflated negative binomial probability function.
Definition: lgamma.hpp:192
Type lgamma(Type x)
Logarithm of gamma function (following R argument convention).
Definition: lgamma.hpp:12
#define VECTORIZE4_ttti(FUN)
Vectorize 4-argument functions.
Definition: Vectorize.hpp:105
Type dnbinom(const Type &x, const Type &size, const Type &prob, int give_log=0)
Negative binomial probability function.Parameterized through size and prob parameters, following R-convention.
Definition: lgamma.hpp:69
Type dnbinom_robust(const Type &x, const Type &log_mu, const Type &log_var_minus_mu, int give_log=0)
Negative binomial probability function.
Definition: lgamma.hpp:103
Type dnbinom2(const Type &x, const Type &mu, const Type &var, int give_log=0)
Negative binomial probability function.Alternative parameterization through mean and variance paramet...
Definition: lgamma.hpp:86
Type dzipois(const Type &x, const Type &lambda, const Type &zip, int give_log=0)
Zero-Inflated Poisson probability function.
Definition: lgamma.hpp:158
Type dzinbinom(const Type &x, const Type &size, const Type &p, const Type &zip, int give_log=0)
Zero-Inflated negative binomial probability function.
Definition: lgamma.hpp:175
Type rnbinom2(Type mu, Type var)
Simulate from a negative binomial distribution.
Definition: lgamma.hpp:218
CppAD::vector< Type > D_lgamma(CppAD::vector< Type > x)
Atomic version of the n&#39;th order derivative of the log gamma function. where the 2 input parameters ...
#define VECTORIZE3_tti(FUN)
Vectorize 3-argument functions.
Definition: Vectorize.hpp:81
License: GPL v2