Stan  1.1
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 #include <stan/agrad.hpp>
5 #include <stan/math/error_handling.hpp>
6 #include <stan/math/special_functions.hpp>
7 #include <stan/meta/traits.hpp>
8 #include <stan/prob/traits.hpp>
9 #include <stan/prob/constants.hpp>
10 #include <stan/prob/internal_math.hpp>
11 
12 namespace stan {
13 
14  namespace prob {
15 
16  // BetaBinomial(n|alpha,beta) [alpha > 0; beta > 0; n >= 0]
17  template <bool propto,
18  typename T_n,
19  typename T_N,
20  typename T_size1,
21  typename T_size2,
22  class Policy>
23  typename return_type<T_size1,T_size2>::type
24  beta_binomial_log(const T_n& n,
25  const T_N& N,
26  const T_size1& alpha,
27  const T_size2& beta,
28  const Policy&) {
29  static const char* function = "stan::prob::beta_binomial_log(%1%)";
30 
31  using stan::math::check_finite;
32  using stan::math::check_nonnegative;
33  using stan::math::check_positive;
34  using stan::math::value_of;
35  using stan::math::check_consistent_sizes;
36  using stan::prob::include_summand;
37 
38  // check if any vectors are zero length
39  if (!(stan::length(n)
40  && stan::length(N)
41  && stan::length(alpha)
42  && stan::length(beta)))
43  return 0.0;
44 
45  typename return_type<T_size1,T_size2>::type logp(0.0);
46  if (!check_nonnegative(function, N, "Population size parameter", &logp, Policy()))
47  return logp;
48  if (!check_finite(function, alpha, "First prior sample size parameter", &logp,
49  Policy()))
50  return logp;
51  if (!check_positive(function, alpha, "First prior sample size parameter",
52  &logp, Policy()))
53  return logp;
54  if (!check_finite(function, beta, "Second prior sample size parameter",
55  &logp, Policy()))
56  return logp;
57  if (!check_positive(function, beta, "Second prior sample size parameter",
58  &logp, Policy()))
59  return logp;
60  if (!(check_consistent_sizes(function,
61  n,N,alpha,beta,
62  "Successes variable","Population size parameter","First prior sample size parameter","Second prior sample size parameter",
63  &logp, Policy())))
64  return logp;
65 
66  // check if no variables are involved and prop-to
67  if (!include_summand<propto,T_size1,T_size2>::value)
68  return 0.0;
69 
70  VectorView<const T_n> n_vec(n);
71  VectorView<const T_N> N_vec(N);
72  VectorView<const T_size1> alpha_vec(alpha);
73  VectorView<const T_size2> beta_vec(beta);
74  size_t size = max_size(n, N, alpha, beta);
75 
76  for (size_t i = 0; i < size; i++) {
77  if (n_vec[i] < 0 || n_vec[i] > N_vec[i])
78  return LOG_ZERO;
79  }
80 
81  using stan::math::lbeta;
82  using stan::math::binomial_coefficient_log;
83 
84  for (size_t i = 0; i < size; i++) {
85  if (include_summand<propto>::value)
86  logp += binomial_coefficient_log(N_vec[i],n_vec[i]);
87  if (include_summand<propto,T_size1,T_size2>::value)
88  logp += lbeta(n_vec[i] + alpha_vec[i], N_vec[i] - n_vec[i] + beta_vec[i])
89  - lbeta(alpha_vec[i],beta_vec[i]);
90  }
91  return logp;
92  }
93 
94  template <bool propto,
95  typename T_n,
96  typename T_N,
97  typename T_size1,
98  typename T_size2>
99  typename return_type<T_size1,T_size2>::type
100  beta_binomial_log(const T_n& n, const T_N& N,
101  const T_size1& alpha, const T_size2& beta) {
102  return beta_binomial_log<propto>(n,N,alpha,beta,
103  stan::math::default_policy());
104  }
105 
106  template <typename T_n,
107  typename T_N,
108  typename T_size1,
109  typename T_size2,
110  class Policy>
111  typename return_type<T_size1,T_size2>::type
112  inline
113  beta_binomial_log(const T_n& n, const T_N& N,
114  const T_size1& alpha, const T_size2& beta,
115  const Policy&) {
116  return beta_binomial_log<false>(n,N,alpha,beta,Policy());
117  }
118 
119  template <typename T_n,
120  typename T_N,
121  typename T_size1,
122  typename T_size2>
123  typename return_type<T_size1,T_size2>::type
124  beta_binomial_log(const T_n& n, const T_N& N,
125  const T_size1& alpha, const T_size2& beta) {
126  return beta_binomial_log<false>(n,N,alpha,beta,
127  stan::math::default_policy());
128  }
129 
130  // Beta-Binomial CDF
131  template <bool propto, typename T_n, typename T_N, typename T_size1, typename T_size2, class Policy>
132  typename return_type<T_size1,T_size2>::type
133  beta_binomial_cdf(const T_n& n, const T_N& N, const T_size1& alpha, const T_size2& beta, const Policy&) {
134 
135  static const char* function = "stan::prob::beta_binomial_cdf(%1%)";
136 
137  using stan::math::check_finite;
138  using stan::math::check_nonnegative;
139  using stan::math::check_positive;
140  using stan::math::value_of;
141  using stan::math::check_consistent_sizes;
142  using stan::prob::include_summand;
143 
144  // Ensure non-zero argument lengths
145  if (!(stan::length(n) && stan::length(N) && stan::length(alpha) && stan::length(beta)))
146  return 0.0;
147 
148  double P(1.0);
149 
150  // Validate arguments
151  if (!check_nonnegative(function, N, "Population size parameter", &P, Policy()))
152  return P;
153 
154  if (!check_finite(function, alpha, "First prior sample size parameter", &P,
155  Policy()))
156  return P;
157 
158  if (!check_positive(function, alpha, "First prior sample size parameter",
159  &P, Policy()))
160  return P;
161 
162  if (!check_finite(function, beta, "Second prior sample size parameter",
163  &P, Policy()))
164  return P;
165 
166  if (!check_positive(function, beta, "Second prior sample size parameter",
167  &P, Policy()))
168  return P;
169 
170  if (!(check_consistent_sizes(function,
171  n, N, alpha, beta,
172  "Successes variable", "Population size parameter",
173  "First prior sample size parameter", "Second prior sample size parameter",
174  &P, Policy())))
175  return P;
176 
177  // Return if everything is constant and only proportionality is required
178  if (!include_summand<propto, T_size1, T_size2>::value)
179  return 0.0;
180 
181  // Wrap arguments in vector views
182  VectorView<const T_n> n_vec(n);
183  VectorView<const T_N> N_vec(N);
184  VectorView<const T_size1> alpha_vec(alpha);
185  VectorView<const T_size2> beta_vec(beta);
186  size_t size = max_size(n, N, alpha, beta);
187 
188  for (size_t i = 0; i < size; i++) {
189 
190  // Return if any of the vectorized trials are inconsistent
191  if (n_vec[i] < 0 || n_vec[i] > N_vec[i])
192  return LOG_ZERO;
193 
194  }
195 
196  // Compute vectorized CDF and gradient
197  using stan::math::value_of;
198  using boost::math::lgamma;
199  using boost::math::digamma;
200 
201  agrad::OperandsAndPartials<T_size1, T_size2> operands_and_partials(alpha, beta);
202 
203  std::fill(operands_and_partials.all_partials,
204  operands_and_partials.all_partials + operands_and_partials.nvaris, 0.0);
205 
206  for (size_t i = 0; i < size; i++) {
207 
208  const double n_dbl = value_of(n_vec[i]);
209  const double N_dbl = value_of(N_vec[i]);
210  const double alpha_dbl = value_of(alpha_vec[i]);
211  const double beta_dbl = value_of(beta_vec[i]);
212 
213  const double mu = alpha_dbl + n_dbl + 1;
214  const double nu = beta_dbl + N_dbl - n_dbl - 1;
215 
216  const double F = stan::math::F32(1, mu, -N_dbl + n_dbl + 1, n_dbl + 2, 1 - nu, 1);
217 
218  double C = lgamma(nu) - lgamma(N_dbl - n_dbl);
219  C += lgamma(mu) - lgamma(n_dbl + 2);
220  C += lgamma(N_dbl + 2) - lgamma(N_dbl + alpha_dbl + beta_dbl);
221  C = std::exp(C);
222 
223  C *= F / boost::math::beta(alpha_dbl, beta_dbl);
224  C /= N_dbl + 1;
225 
226  const double Pi = 1 - C;
227 
228  P *= Pi;
229 
230  double dF[6];
231  double digammaOne = 0;
232  double digammaTwo = 0;
233 
234  if ( (!is_constant_struct<T_size1>::value) || (!is_constant_struct<T_size2>::value) ) {
235 
236  digammaOne = digamma(mu + nu);
237  digammaTwo = digamma(alpha_dbl + beta_dbl);
238 
239  stan::math::gradF32(dF, 1, mu, -N_dbl + n_dbl + 1, n_dbl + 2, 1 - nu, 1);
240 
241  }
242 
243  if (!is_constant_struct<T_size1>::value) {
244 
245  const double g = - C * (digamma(mu) - digammaOne + dF[1] / F - digamma(alpha_dbl) + digammaTwo);
246 
247  operands_and_partials.d_x1[n]
248  += g / Pi;
249 
250  }
251 
252  if (!is_constant_struct<T_size2>::value) {
253 
254  const double g = - C * (digamma(nu) - digammaOne - dF[4] / F - digamma(beta_dbl) + digammaTwo);
255 
256  operands_and_partials.d_x2[n]
257  += g / Pi;
258 
259  }
260 
261  }
262 
263  for (size_t i = 0; i < size; i++) {
264 
265  if (!is_constant_struct<T_size1>::value)
266  operands_and_partials.d_x1[i] *= P;
267 
268  if (!is_constant_struct<T_size2>::value)
269  operands_and_partials.d_x2[i] *= P;
270 
271  }
272 
273  return P;
274 
275  }
276 
277  template <bool propto, typename T_n, typename T_N, typename T_size1, typename T_size2>
278  inline typename return_type<T_size1,T_size2>::type
279  beta_binomial_cdf(const T_n& n, const T_N& N, const T_size1& alpha, const T_size2& beta) {
280  return beta_binomial_cdf<propto>(n, N, alpha, beta, stan::math::default_policy());
281  }
282 
283  template <typename T_n, typename T_N, typename T_size1, typename T_size2, class Policy>
284  inline typename return_type<T_size1,T_size2>::type
285  beta_binomial_cdf(const T_n& n, const T_N& N, const T_size1& alpha, const T_size2& beta, const Policy&) {
286  return beta_binomial_cdf<false>(n, N, alpha, beta, Policy());
287  }
288 
289  template <typename T_n, typename T_N, typename T_size1, typename T_size2>
290  typename return_type<T_size1,T_size2>::type
291  beta_binomial_cdf(const T_n& n, const T_N& N, const T_size1& alpha, const T_size2& beta) {
292  return beta_binomial_cdf<false>(n, N, alpha, beta, stan::math::default_policy());
293  }
294 
295 
296  }
297 }
298 #endif
     [ Stan Home Page ] © 2011–2013, Stan Development Team.