Stan  1.1
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
neg_binomial.hpp
Go to the documentation of this file.
1 #ifndef __STAN__PROB__DISTRIBUTIONS__UNIVARIATE__DISCRETE__NEG_BINOMIAL_HPP__
2 #define __STAN__PROB__DISTRIBUTIONS__UNIVARIATE__DISCRETE__NEG_BINOMIAL_HPP__
3 
4 #include <boost/math/special_functions/digamma.hpp>
5 #include <stan/agrad.hpp>
6 #include <stan/math/error_handling.hpp>
7 #include <stan/math/special_functions.hpp>
8 #include <stan/meta/traits.hpp>
9 #include <stan/prob/traits.hpp>
10 #include <stan/prob/constants.hpp>
11 #include <stan/prob/internal_math.hpp>
12 
13 namespace stan {
14 
15  namespace prob {
16 
17  // NegBinomial(n|alpha,beta) [alpha > 0; beta > 0; n >= 0]
18  template <bool propto,
19  typename T_n,
20  typename T_shape, typename T_inv_scale,
21  class Policy>
22  typename return_type<T_shape, T_inv_scale>::type
23  neg_binomial_log(const T_n& n,
24  const T_shape& alpha,
25  const T_inv_scale& beta,
26  const Policy&) {
27 
28  static const char* function = "stan::prob::neg_binomial_log(%1%)";
29 
30  using stan::math::check_finite;
31  using stan::math::check_nonnegative;
32  using stan::math::check_positive;
33  using stan::math::value_of;
34  using stan::math::check_consistent_sizes;
35  using stan::prob::include_summand;
36 
37  // check if any vectors are zero length
38  if (!(stan::length(n)
39  && stan::length(alpha)
40  && stan::length(beta)))
41  return 0.0;
42 
43  double logp(0.0);
44  if (!check_nonnegative(function, n, "Failures variable", &logp, Policy()))
45  return logp;
46  if (!check_finite(function, alpha, "Shape parameter", &logp, Policy()))
47  return logp;
48  if (!check_positive(function, alpha, "Shape parameter", &logp, Policy()))
49  return logp;
50  if (!check_finite(function, beta, "Inverse scale parameter",
51  &logp, Policy()))
52  return logp;
53  if (!check_positive(function, beta, "Inverse scale parameter",
54  &logp, Policy()))
55  return logp;
56  if (!(check_consistent_sizes(function,
57  n,alpha,beta,
58  "Failures variable",
59  "Shape parameter","Inverse scale parameter",
60  &logp, Policy())))
61  return logp;
62 
63  // check if no variables are involved and prop-to
64  if (!include_summand<propto,T_shape,T_inv_scale>::value)
65  return 0.0;
66 
67  using stan::math::multiply_log;
68  using stan::math::binomial_coefficient_log;
69  using boost::math::digamma;
70 
71  // set up template expressions wrapping scalars into vector views
72  VectorView<const T_n> n_vec(n);
73  VectorView<const T_shape> alpha_vec(alpha);
74  VectorView<const T_inv_scale> beta_vec(beta);
75  size_t size = max_size(n, alpha, beta);
76 
77  agrad::OperandsAndPartials<T_shape,T_inv_scale> operands_and_partials(alpha,beta);
78 
79  size_t len_ab = max_size(alpha,beta);
80  DoubleVectorView<true,(is_vector<T_shape>::value || is_vector<T_inv_scale>::value)>
81  lambda(len_ab);
82  for (size_t i = 0; i < len_ab; ++i)
83  lambda[i] = value_of(alpha_vec[i]) / value_of(beta_vec[i]);
84 
85  DoubleVectorView<true,is_vector<T_inv_scale>::value>
86  log1p_beta(length(beta));
87  for (size_t i = 0; i < length(beta); ++i)
88  log1p_beta[i] = log1p(value_of(beta_vec[i]));
89 
90  DoubleVectorView<true,is_vector<T_inv_scale>::value>
91  log_beta_m_log1p_beta(length(beta));
92  for (size_t i = 0; i < length(beta); ++i)
93  log_beta_m_log1p_beta[i] = log(value_of(beta_vec[i])) - log1p_beta[i];
94 
95  DoubleVectorView<true,(is_vector<T_inv_scale>::value || is_vector<T_shape>::value)>
96  alpha_times_log_beta_over_1p_beta(len_ab);
97  for (size_t i = 0; i < len_ab; ++i)
98  alpha_times_log_beta_over_1p_beta[i]
99  = value_of(alpha_vec[i])
100  * log(value_of(beta_vec[i])
101  / (1.0 + value_of(beta_vec[i])));
102 
103  DoubleVectorView<!is_constant_struct<T_shape>::value, is_vector<T_shape>::value>
104  digamma_alpha(length(alpha));
105  if (!is_constant_struct<T_shape>::value)
106  for (size_t i = 0; i < length(alpha); ++i)
107  digamma_alpha[i] = digamma(value_of(alpha_vec[i]));
108 
109  DoubleVectorView<!is_constant_struct<T_shape>::value, is_vector<T_inv_scale>::value>
110  log_beta(length(beta));
111  if (!is_constant_struct<T_shape>::value)
112  for (size_t i = 0; i < length(beta); ++i)
113  log_beta[i] = log(value_of(beta_vec[i]));
114 
115  DoubleVectorView<!is_constant_struct<T_inv_scale>::value,
116  (is_vector<T_shape>::value || is_vector<T_inv_scale>::value)>
117  lambda_m_alpha_over_1p_beta(len_ab);
118  if (!is_constant_struct<T_inv_scale>::value)
119  for (size_t i = 0; i < len_ab; ++i)
120  lambda_m_alpha_over_1p_beta[i] =
121  lambda[i]
122  - ( value_of(alpha_vec[i])
123  / (1.0 + value_of(beta_vec[i])) );
124 
125  for (size_t i = 0; i < size; i++) {
126  if (alpha_vec[i] > 1e10) { // reduces numerically to Poisson
127  if (include_summand<propto>::value)
128  logp -= lgamma(n_vec[i] + 1.0);
129  if (include_summand<propto,T_shape,T_inv_scale>::value)
130  logp += multiply_log(n_vec[i], lambda[i]) - lambda[i];
131 
132  if (!is_constant_struct<T_shape>::value)
133  operands_and_partials.d_x1[i]
134  += n_vec[i] / value_of(alpha_vec[i])
135  - 1.0 / value_of(beta_vec[i]);
136  if (!is_constant_struct<T_inv_scale>::value)
137  operands_and_partials.d_x2[i]
138  += (lambda[i] - n_vec[i]) / value_of(beta_vec[i]) ;
139  } else { // standard density definition
140  if (include_summand<propto,T_shape>::value)
141  if (n_vec[i] != 0)
142  logp += binomial_coefficient_log<double>(n_vec[i]
143  + value_of(alpha_vec[i]) - 1.0,
144  n_vec[i]);
145  if (include_summand<propto,T_shape,T_inv_scale>::value)
146  logp +=
147  alpha_times_log_beta_over_1p_beta[i]
148  - n_vec[i] * log1p_beta[i];
149 
150  if (!is_constant_struct<T_shape>::value)
151  operands_and_partials.d_x1[i]
152  += digamma(value_of(alpha_vec[i]) + n_vec[i])
153  - digamma_alpha[i]
154  + log_beta_m_log1p_beta[i];
155  if (!is_constant_struct<T_inv_scale>::value)
156  operands_and_partials.d_x2[i]
157  += lambda_m_alpha_over_1p_beta[i]
158  - n_vec[i] / (value_of(beta_vec[i]) + 1.0);
159  }
160  }
161  return operands_and_partials.to_var(logp);
162  }
163 
164  template <bool propto,
165  typename T_n,
166  typename T_shape, typename T_inv_scale>
167  inline
168  typename return_type<T_shape, T_inv_scale>::type
169  neg_binomial_log(const T_n& n,
170  const T_shape& alpha,
171  const T_inv_scale& beta) {
172  return neg_binomial_log<propto>(n,alpha,beta,
173  stan::math::default_policy());
174  }
175 
176  template <typename T_n,
177  typename T_shape, typename T_inv_scale,
178  class Policy>
179  inline
180  typename return_type<T_shape, T_inv_scale>::type
181  neg_binomial_log(const T_n& n,
182  const T_shape& alpha,
183  const T_inv_scale& beta,
184  const Policy&) {
185  return neg_binomial_log<false>(n,alpha,beta,Policy());
186  }
187 
188  template <typename T_n,
189  typename T_shape, typename T_inv_scale>
190  inline
191  typename return_type<T_shape, T_inv_scale>::type
192  neg_binomial_log(const T_n& n,
193  const T_shape& alpha,
194  const T_inv_scale& beta) {
195  return neg_binomial_log<false>(n,alpha,beta,
196  stan::math::default_policy());
197  }
198 
199  // Negative Binomial CDF
200  template <bool propto, typename T_n, typename T_shape, typename T_inv_scale,
201  class Policy>
202  typename return_type<T_shape, T_inv_scale>::type
203  neg_binomial_cdf(const T_n& n, const T_shape& alpha, const T_inv_scale& beta,
204  const Policy&) {
205 
206  static const char* function = "stan::prob::neg_binomial_cdf(%1%)";
207 
208  using stan::math::check_finite;
209  using stan::math::check_nonnegative;
210  using stan::math::check_positive;
211  using stan::math::check_consistent_sizes;
212  using stan::prob::include_summand;
213 
214  // Ensure non-zero arugment lengths
215  if (!(stan::length(n) && stan::length(alpha) && stan::length(beta)))
216  return 0.0;
217 
218  double P(1.0);
219 
220  // Validate arguments
221  if (!check_nonnegative(function, n, "Failures variable", &P, Policy()))
222  return P;
223 
224  if (!check_finite(function, alpha, "Shape parameter", &P, Policy()))
225  return P;
226 
227  if (!check_positive(function, alpha, "Shape parameter", &P, Policy()))
228  return P;
229 
230  if (!check_finite(function, beta, "Inverse scale parameter",
231  &P, Policy()))
232  return P;
233 
234  if (!check_positive(function, beta, "Inverse scale parameter",
235  &P, Policy()))
236  return P;
237 
238  if (!(check_consistent_sizes(function,
239  n, alpha, beta,
240  "Failures variable",
241  "Shape parameter","Inverse scale parameter",
242  &P, Policy())))
243  return P;
244 
245  // Return if everything is constant and only proportionality is required
246  if (!include_summand<propto,T_shape,T_inv_scale>::value)
247  return 0.0;
248 
249  // Wrap arguments in vector views
250  VectorView<const T_n> n_vec(n);
251  VectorView<const T_shape> alpha_vec(alpha);
252  VectorView<const T_inv_scale> beta_vec(beta);
253  size_t size = max_size(n, alpha, beta);
254 
255  // Compute vectorized CDF and gradient
256  using stan::math::value_of;
257  using boost::math::ibetac;
258  using boost::math::ibeta_derivative;
259 
260  using boost::math::digamma;
261 
262  // Cache a few expensive function calls if alpha is a parameter
263  DoubleVectorView<!is_constant_struct<T_shape>::value,
264  is_vector<T_shape>::value> digammaN_vec(stan::length(alpha));
265  DoubleVectorView<!is_constant_struct<T_shape>::value,
266  is_vector<T_shape>::value> digammaAlpha_vec(stan::length(alpha));
267  DoubleVectorView<!is_constant_struct<T_shape>::value,
268  is_vector<T_shape>::value> digammaSum_vec(stan::length(alpha));
269  DoubleVectorView<!is_constant_struct<T_shape>::value,
270  is_vector<T_shape>::value> betaFunc_vec(stan::length(alpha));
271 
272  if (!is_constant_struct<T_shape>::value) {
273  for (size_t i = 0; i < stan::length(alpha); i++) {
274  const double n_dbl = value_of(n_vec[i]);
275  const double alpha_dbl = value_of(alpha_vec[i]);
276 
277  digammaN_vec[i] = digamma(n_dbl);
278  digammaAlpha_vec[i] = digamma(alpha_dbl + 1);
279  digammaSum_vec[i] = digamma(n_dbl + alpha_dbl + 1);
280  betaFunc_vec[i] = boost::math::beta(n_dbl, alpha_dbl + 1);
281  }
282  }
283 
284  agrad::OperandsAndPartials<T_shape, T_inv_scale>
285  operands_and_partials(alpha, beta);
286 
287  std::fill(operands_and_partials.all_partials,
288  operands_and_partials.all_partials + operands_and_partials.nvaris,
289  0.0);
290 
291  for (size_t i = 0; i < size; i++) {
292  const double n_dbl = value_of(n_vec[i]);
293  const double alpha_dbl = value_of(alpha_vec[i]);
294  const double beta_dbl = value_of(beta_vec[i]);
295 
296  const double p_dbl = beta_dbl / (1.0 + beta_dbl);
297  const double d_dbl = 1.0 / ( (1.0 + beta_dbl) * (1.0 + beta_dbl) );
298 
299  const double Pi = ibetac(alpha_dbl + 1, n_dbl, p_dbl);
300 
301  P *= Pi;
302 
303  if (!is_constant_struct<T_shape>::value) {
304  double g1 = 0;
305  double g2 = 0;
306 
307  stan::math::gradRegIncBeta(g1, g2, alpha_dbl + 1, n_dbl, p_dbl,
308  digammaAlpha_vec[i], digammaN_vec[i],
309  digammaSum_vec[i], betaFunc_vec[i]);
310 
311  operands_and_partials.d_x1[n]
312  += - g1 / Pi;
313  }
314 
315  if (!is_constant_struct<T_inv_scale>::value)
316  operands_and_partials.d_x2[i]
317  += - d_dbl * ibeta_derivative(alpha_dbl + 1, n_dbl, p_dbl) / Pi;
318 
319  }
320 
321  for (size_t i = 0; i < size; i++) {
322  if (!is_constant_struct<T_shape>::value)
323  operands_and_partials.d_x1[i] *= P;
324  if (!is_constant_struct<T_inv_scale>::value)
325  operands_and_partials.d_x2[i] *= P;
326  }
327 
328  return P;
329  }
330 
331  template <bool propto, typename T_n, typename T_shape, typename T_inv_scale>
332  inline typename return_type<T_shape, T_inv_scale>::type
333  neg_binomial_cdf(const T_n& n, const T_shape& alpha, const T_inv_scale& beta) {
334  return neg_binomial_cdf<propto>(n, alpha, beta, stan::math::default_policy());
335  }
336 
337  template <typename T_n, typename T_shape, typename T_inv_scale, class Policy>
338  inline typename return_type<T_shape, T_inv_scale>::type
339  neg_binomial_cdf(const T_n& n, const T_shape& alpha, const T_inv_scale& beta,
340  const Policy&) {
341  return neg_binomial_cdf<false>(n, alpha, beta, Policy());
342  }
343 
344  template <typename T_n, typename T_shape, typename T_inv_scale>
345  inline typename return_type<T_shape, T_inv_scale>::type
346  neg_binomial_cdf(const T_n& n, const T_shape& alpha, const T_inv_scale& beta) {
347  return neg_binomial_cdf<false>(n, alpha, beta, stan::math::default_policy());
348  }
349 
350 
351  }
352 }
353 #endif
     [ Stan Home Page ] © 2011–2013, Stan Development Team.