Stan  1.3
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
beta_binomial.hpp
Go to the documentation of this file.
1 #ifndef __STAN__PROB__DISTRIBUTIONS__UNIVARIATE__DISCRETE__BETA_BINOMIAL_HPP__
2 #define __STAN__PROB__DISTRIBUTIONS__UNIVARIATE__DISCRETE__BETA_BINOMIAL_HPP__
3 
4 
5 #include <stan/prob/distributions/univariate/discrete/binomial.hpp>
6 #include <stan/prob/distributions/univariate/continuous/beta.hpp>
7 
8 #include <stan/agrad.hpp>
9 #include <stan/math/error_handling.hpp>
10 #include <stan/math/functions/lbeta.hpp>
11 #include <stan/math/functions/value_of.hpp>
12 #include <stan/meta/traits.hpp>
13 #include <stan/prob/traits.hpp>
14 #include <stan/prob/constants.hpp>
15 #include <stan/prob/internal_math.hpp>
16 
17 #include <stan/math/functions/binomial_coefficient_log.hpp>
18 
19 namespace stan {
20 
21  namespace prob {
22 
23  // BetaBinomial(n|alpha,beta) [alpha > 0; beta > 0; n >= 0]
24  template <bool propto,
25  typename T_n,
26  typename T_N,
27  typename T_size1,
28  typename T_size2,
29  class Policy>
30  typename return_type<T_size1,T_size2>::type
31  beta_binomial_log(const T_n& n,
32  const T_N& N,
33  const T_size1& alpha,
34  const T_size2& beta,
35  const Policy&) {
36  static const char* function = "stan::prob::beta_binomial_log(%1%)";
37 
38  using stan::math::check_finite;
39  using stan::math::check_nonnegative;
40  using stan::math::check_positive;
41  using stan::math::value_of;
42  using stan::math::check_consistent_sizes;
43  using stan::prob::include_summand;
44 
45  // check if any vectors are zero length
46  if (!(stan::length(n)
47  && stan::length(N)
48  && stan::length(alpha)
49  && stan::length(beta)))
50  return 0.0;
51 
52  double logp(0.0);
53  if (!check_nonnegative(function, N, "Population size parameter",
54  &logp, Policy()))
55  return logp;
56  if (!check_finite(function, alpha, "First prior sample size parameter",
57  &logp, Policy()))
58  return logp;
59  if (!check_positive(function, alpha, "First prior sample size parameter",
60  &logp, Policy()))
61  return logp;
62  if (!check_finite(function, beta, "Second prior sample size parameter",
63  &logp, Policy()))
64  return logp;
65  if (!check_positive(function, beta, "Second prior sample size parameter",
66  &logp, Policy()))
67  return logp;
68  if (!(check_consistent_sizes(function,
69  n,N,alpha,beta,
70  "Successes variable",
71  "Population size parameter",
72  "First prior sample size parameter",
73  "Second prior sample size parameter",
74  &logp, Policy())))
75  return logp;
76 
77  // check if no variables are involved and prop-to
78  if (!include_summand<propto,T_size1,T_size2>::value)
79  return 0.0;
80 
81  VectorView<const T_n> n_vec(n);
82  VectorView<const T_N> N_vec(N);
83  VectorView<const T_size1> alpha_vec(alpha);
84  VectorView<const T_size2> beta_vec(beta);
85  size_t size = max_size(n, N, alpha, beta);
86 
87  for (size_t i = 0; i < size; i++) {
88  if (n_vec[i] < 0 || n_vec[i] > N_vec[i])
89  return LOG_ZERO;
90  }
91 
92  using stan::math::lbeta;
93  using stan::math::binomial_coefficient_log;
94  using boost::math::digamma;
95 
96  DoubleVectorView<include_summand<propto>::value,
97  is_vector<T_n>::value || is_vector<T_N>::value>
98  normalizing_constant(max_size(N,n));
99  for (size_t i = 0; i < max_size(N,n); i++)
100  if (include_summand<propto>::value)
101  normalizing_constant[i]
102  = binomial_coefficient_log(N_vec[i],n_vec[i]);
103 
104  DoubleVectorView<include_summand<propto,T_size1,T_size2>::value,
105  is_vector<T_n>::value || is_vector<T_N>::value
106  || is_vector<T_size1>::value || is_vector<T_size2>::value>
107  lbeta_numerator(size);
108  for (size_t i = 0; i < size; i++)
109  if (include_summand<propto,T_size1,T_size2>::value)
110  lbeta_numerator[i] = lbeta(n_vec[i] + value_of(alpha_vec[i]),
111  N_vec[i] - n_vec[i]
112  + value_of(beta_vec[i]));
113  DoubleVectorView<include_summand<propto,T_size1,T_size2>::value,
114  is_vector<T_size1>::value || is_vector<T_size2>::value>
115  lbeta_denominator(max_size(alpha,beta));
116 for (size_t i = 0; i < max_size(alpha,beta); i++)
117  if (include_summand<propto,T_size1,T_size2>::value)
118  lbeta_denominator[i] = lbeta(value_of(alpha_vec[i]),
119  value_of(beta_vec[i]));
120 
121 DoubleVectorView<!is_constant_struct<T_size1>::value,
122  is_vector<T_n>::value || is_vector<T_size1>::value>
123  digamma_n_plus_alpha(max_size(n,alpha));
124  for (size_t i = 0; i < max_size(n,alpha); i++)
125  if (!is_constant_struct<T_size1>::value)
126  digamma_n_plus_alpha[i]
127  = digamma(n_vec[i] + value_of(alpha_vec[i]));
128 
129  DoubleVectorView<!is_constant_struct<T_size1>::value
130  || !is_constant_struct<T_size2>::value,
131  is_vector<T_N>::value
132  || is_vector<T_size1>::value
133  || is_vector<T_size1>::value>
134  digamma_N_plus_alpha_plus_beta(max_size(N,alpha,beta));
135  for (size_t i = 0; i < max_size(N,alpha,beta); i++)
136  if (!is_constant_struct<T_size1>::value
137  || !is_constant_struct<T_size2>::value)
138  digamma_N_plus_alpha_plus_beta[i]
139  = digamma(N_vec[i] + value_of(alpha_vec[i]) + value_of(beta_vec[i]));
140 
141  DoubleVectorView<!is_constant_struct<T_size1>::value
142  || !is_constant_struct<T_size2>::value,
143  is_vector<T_size1>::value
144  || is_vector<T_size1>::value>
145  digamma_alpha_plus_beta(max_size(alpha,beta));
146  for (size_t i = 0; i < max_size(alpha,beta); i++)
147  if (!is_constant_struct<T_size1>::value
148  || !is_constant_struct<T_size2>::value)
149  digamma_alpha_plus_beta[i]
150  = digamma(value_of(alpha_vec[i]) + value_of(beta_vec[i]));
151 
152  DoubleVectorView<!is_constant_struct<T_size1>::value,
153  is_vector<T_size1>::value>
154  digamma_alpha(length(alpha));
155 for (size_t i = 0; i < length(alpha); i++)
156  if (!is_constant_struct<T_size1>::value)
157  digamma_alpha[i] = digamma(value_of(alpha_vec[i]));
158 
159 DoubleVectorView<!is_constant_struct<T_size2>::value,
160  is_vector<T_size2>::value>
161  digamma_beta(length(beta));
162  for (size_t i = 0; i < length(beta); i++)
163  if (!is_constant_struct<T_size2>::value)
164  digamma_beta[i] = digamma(value_of(beta_vec[i]));
165 
166  agrad::OperandsAndPartials<T_n,T_N,T_size1,T_size2>
167  operands_and_partials(n,N,alpha,beta);
168  for (size_t i = 0; i < size; i++) {
169  if (include_summand<propto>::value)
170  logp += normalizing_constant[i];
171  if (include_summand<propto,T_size1,T_size2>::value)
172  logp += lbeta_numerator[i]
173  - lbeta_denominator[i];
174 
175  if (!is_constant_struct<T_size1>::value)
176  operands_and_partials.d_x3[i]
177  += digamma_n_plus_alpha[i]
178  - digamma_N_plus_alpha_plus_beta[i]
179  + digamma_alpha_plus_beta[i]
180  - digamma_alpha[i];
181  if (!is_constant_struct<T_size2>::value)
182  operands_and_partials.d_x4[i]
183  += digamma(value_of(N_vec[i]-n_vec[i]+beta_vec[i]))
184  - digamma_N_plus_alpha_plus_beta[i]
185  + digamma_alpha_plus_beta[i]
186  - digamma_beta[i];
187  }
188  return operands_and_partials.to_var(logp);
189  }
190 
191  template <bool propto,
192  typename T_n,
193  typename T_N,
194  typename T_size1,
195  typename T_size2>
196  typename return_type<T_size1,T_size2>::type
197  beta_binomial_log(const T_n& n, const T_N& N,
198  const T_size1& alpha, const T_size2& beta) {
199  return beta_binomial_log<propto>(n,N,alpha,beta,
200  stan::math::default_policy());
201  }
202 
203  template <typename T_n,
204  typename T_N,
205  typename T_size1,
206  typename T_size2,
207  class Policy>
208  typename return_type<T_size1,T_size2>::type
209  inline
210  beta_binomial_log(const T_n& n, const T_N& N,
211  const T_size1& alpha, const T_size2& beta,
212  const Policy&) {
213  return beta_binomial_log<false>(n,N,alpha,beta,Policy());
214  }
215 
216  template <typename T_n,
217  typename T_N,
218  typename T_size1,
219  typename T_size2>
220  typename return_type<T_size1,T_size2>::type
221  beta_binomial_log(const T_n& n, const T_N& N,
222  const T_size1& alpha, const T_size2& beta) {
223  return beta_binomial_log<false>(n,N,alpha,beta,
224  stan::math::default_policy());
225  }
226 
227  // Beta-Binomial CDF
228  template <bool propto, typename T_n, typename T_N, typename T_size1,
229  typename T_size2, class Policy>
230  typename return_type<T_size1,T_size2>::type
231  beta_binomial_cdf(const T_n& n, const T_N& N, const T_size1& alpha,
232  const T_size2& beta, const Policy&) {
233 
234  static const char* function = "stan::prob::beta_binomial_cdf(%1%)";
235 
236  using stan::math::check_finite;
237  using stan::math::check_nonnegative;
238  using stan::math::check_positive;
239  using stan::math::value_of;
240  using stan::math::check_consistent_sizes;
241  using stan::prob::include_summand;
242 
243  // Ensure non-zero argument lengths
244  if (!(stan::length(n) && stan::length(N) && stan::length(alpha)
245  && stan::length(beta)))
246  return 1.0;
247 
248  double P(1.0);
249 
250  // Validate arguments
251  if (!check_nonnegative(function, N, "Population size parameter",
252  &P, Policy()))
253  return P;
254 
255  if (!check_finite(function, alpha, "First prior sample size parameter",
256  &P, Policy()))
257  return P;
258 
259  if (!check_positive(function, alpha, "First prior sample size parameter",
260  &P, Policy()))
261  return P;
262 
263  if (!check_finite(function, beta, "Second prior sample size parameter",
264  &P, Policy()))
265  return P;
266 
267  if (!check_positive(function, beta, "Second prior sample size parameter",
268  &P, Policy()))
269  return P;
270 
271  if (!(check_consistent_sizes(function,
272  n, N, alpha, beta,
273  "Successes variable",
274  "Population size parameter",
275  "First prior sample size parameter",
276  "Second prior sample size parameter",
277  &P, Policy())))
278  return P;
279 
280  // Return if everything is constant and only proportionality is required
281  if (!include_summand<propto, T_size1, T_size2>::value)
282  return 1.0;
283 
284  // Wrap arguments in vector views
285  VectorView<const T_n> n_vec(n);
286  VectorView<const T_N> N_vec(N);
287  VectorView<const T_size1> alpha_vec(alpha);
288  VectorView<const T_size2> beta_vec(beta);
289  size_t size = max_size(n, N, alpha, beta);
290 
291  // Compute vectorized CDF and gradient
292  using boost::math::lgamma;
293  using boost::math::digamma;
294 
295  agrad::OperandsAndPartials<T_size1, T_size2>
296  operands_and_partials(alpha, beta);
297 
298  std::fill(operands_and_partials.all_partials,
299  operands_and_partials.all_partials
300  + operands_and_partials.nvaris,
301  0.0);
302 
303  // Explicit return for extreme values
304  // The gradients are technically ill-defined, but treated as zero
305  for (size_t i = 0; i < stan::length(n); i++) {
306  if (value_of(n_vec[i]) <= 0)
307  return operands_and_partials.to_var(0.0);
308  }
309 
310  for (size_t i = 0; i < size; i++) {
311  // Explicit results for extreme values
312  // The gradients are technically ill-defined, but treated as zero
313  if (value_of(n_vec[i]) >= value_of(N_vec[i])) {
314  continue;
315  }
316 
317  const double n_dbl = value_of(n_vec[i]);
318  const double N_dbl = value_of(N_vec[i]);
319  const double alpha_dbl = value_of(alpha_vec[i]);
320  const double beta_dbl = value_of(beta_vec[i]);
321 
322  const double mu = alpha_dbl + n_dbl + 1;
323  const double nu = beta_dbl + N_dbl - n_dbl - 1;
324 
325  const double F = stan::math::F32(1, mu, -N_dbl + n_dbl + 1, n_dbl + 2,
326  1 - nu, 1);
327 
328  double C = lgamma(nu) - lgamma(N_dbl - n_dbl);
329  C += lgamma(mu) - lgamma(n_dbl + 2);
330  C += lgamma(N_dbl + 2) - lgamma(N_dbl + alpha_dbl + beta_dbl);
331  C = std::exp(C);
332 
333  C *= F / boost::math::beta(alpha_dbl, beta_dbl);
334  C /= N_dbl + 1;
335 
336  const double Pi = 1 - C;
337 
338  P *= Pi;
339 
340  double dF[6];
341  double digammaOne = 0;
342  double digammaTwo = 0;
343 
344  if ( (!is_constant_struct<T_size1>::value)
345  || (!is_constant_struct<T_size2>::value) ) {
346 
347  digammaOne = digamma(mu + nu);
348  digammaTwo = digamma(alpha_dbl + beta_dbl);
349 
350  stan::math::gradF32(dF, 1, mu, -N_dbl + n_dbl + 1, n_dbl + 2,
351  1 - nu, 1);
352  }
353 
354  if (!is_constant_struct<T_size1>::value) {
355 
356  const double g
357  = - C * (digamma(mu) - digammaOne + dF[1] / F
358  - digamma(alpha_dbl) + digammaTwo);
359 
360  operands_and_partials.d_x1[i]
361  += g / Pi;
362  }
363 
364  if (!is_constant_struct<T_size2>::value) {
365 
366  const double g
367  = - C * (digamma(nu) - digammaOne - dF[4] / F - digamma(beta_dbl)
368  + digammaTwo);
369 
370  operands_and_partials.d_x2[i]
371  += g / Pi;
372  }
373  }
374 
375  if (!is_constant_struct<T_size1>::value)
376  for(size_t i = 0; i < stan::length(alpha); ++i)
377  operands_and_partials.d_x1[i] *= P;
378 
379 
380  if (!is_constant_struct<T_size2>::value)
381  for(size_t i = 0; i < stan::length(beta); ++i)
382  operands_and_partials.d_x2[i] *= P;
383 
384  return operands_and_partials.to_var(P);
385  }
386 
387  template <bool propto, typename T_n, typename T_N, typename T_size1,
388  typename T_size2>
389  inline typename return_type<T_size1,T_size2>::type
390  beta_binomial_cdf(const T_n& n, const T_N& N, const T_size1& alpha,
391  const T_size2& beta) {
392  return beta_binomial_cdf<propto>(n, N, alpha, beta,
393  stan::math::default_policy());
394  }
395 
396  template <typename T_n, typename T_N, typename T_size1, typename T_size2,
397  class Policy>
398  inline typename return_type<T_size1,T_size2>::type
399  beta_binomial_cdf(const T_n& n, const T_N& N, const T_size1& alpha,
400  const T_size2& beta, const Policy&) {
401  return beta_binomial_cdf<false>(n, N, alpha, beta, Policy());
402  }
403 
404  template <typename T_n, typename T_N, typename T_size1, typename T_size2>
405  typename return_type<T_size1,T_size2>::type
406  beta_binomial_cdf(const T_n& n, const T_N& N, const T_size1& alpha,
407  const T_size2& beta) {
408  return beta_binomial_cdf<false>(n, N, alpha, beta,
409  stan::math::default_policy());
410  }
411 
412  template <class RNG>
413  inline int
414  beta_binomial_rng(const int N,
415  const double alpha,
416  const double beta,
417  RNG& rng) {
418  double a = stan::prob::beta_rng(alpha, beta, rng);
419  while(a > 1 || a < 0)
420  a = stan::prob::beta_rng(alpha, beta, rng);
421  return stan::prob::binomial_rng(N, a, rng);
422  }
423  }
424 }
425 #endif
     [ Stan Home Page ] © 2011–2013, Stan Development Team.