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