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