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