Stan  1.1
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
bernoulli.hpp
Go to the documentation of this file.
1 #ifndef __STAN__PROB__DISTRIBUTIONS__UNIVARIATE__DISCRETE__BERNOULLI_HPP__
2 #define __STAN__PROB__DISTRIBUTIONS__UNIVARIATE__DISCRETE__BERNOULLI_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  // Bernoulli(n|theta) [0 <= n <= 1; 0 <= theta <= 1]
16  // FIXME: documentation
17  template <bool propto,
18  typename T_n, typename T_prob,
19  class Policy>
20  typename return_type<T_prob>::type
21  bernoulli_log(const T_n& n,
22  const T_prob& theta,
23  const Policy&) {
24  static const char* function = "stan::prob::bernoulli_log(%1%)";
25 
26  using stan::math::check_finite;
27  using stan::math::check_bounded;
28  using stan::math::log1m;
29  using stan::math::value_of;
30  using stan::math::check_consistent_sizes;
31  using stan::prob::include_summand;
32 
33  // check if any vectors are zero length
34  if (!(stan::length(n)
35  && stan::length(theta)))
36  return 0.0;
37 
38  // set up return value accumulator
39  double logp(0.0);
40 
41  // validate args (here done over var, which should be OK)
42  if (!check_bounded(function, n, 0, 1, "n", &logp, Policy()))
43  return logp;
44  if (!check_finite(function, theta, "Probability parameter", &logp, Policy()))
45  return logp;
46  if (!check_bounded(function, theta, 0.0, 1.0,
47  "Probability parameter", &logp, Policy()))
48  return logp;
49  if (!(check_consistent_sizes(function,
50  n,theta,
51  "Random variable","Probability parameter",
52  &logp, Policy())))
53  return logp;
54 
55  // check if no variables are involved and prop-to
56  if (!include_summand<propto,T_prob>::value)
57  return 0.0;
58 
59  // set up template expressions wrapping scalars into vector views
60  VectorView<const T_n> n_vec(n);
61  VectorView<const T_prob> theta_vec(theta);
62  size_t N = max_size(n, theta);
63  agrad::OperandsAndPartials<T_prob> operands_and_partials(theta);
64 
65  if (length(theta) == 1) {
66  size_t sum = 0;
67  for (size_t n = 0; n < N; n++) {
68  sum += value_of(n_vec[n]);
69  }
70  const double theta_dbl = value_of(theta_vec[0]);
71  // avoid nans when sum == N or sum == 0
72  if (sum == N) {
73  logp += N * log(theta_dbl);
74  operands_and_partials.d_x1[0] += N / theta_dbl;
75  } else if (sum == 0) {
76  logp += N * log1m(theta_dbl);
77  operands_and_partials.d_x1[0] += N / (theta_dbl - 1);
78  } else {
79  const double log_theta = log(theta_dbl);
80  const double log1m_theta = log1m(theta_dbl);
81  if (include_summand<propto,T_prob>::value) {
82  logp += sum * log_theta;
83  logp += (N - sum) * log1m_theta;
84 
85  operands_and_partials.d_x1[0] += sum / theta_dbl;
86  operands_and_partials.d_x1[0] += (N - sum) / (theta_dbl - 1);
87  }
88  }
89  } else {
90  for (size_t n = 0; n < N; n++) {
91  // pull out values of arguments
92  const int n_int = value_of(n_vec[n]);
93  const double theta_dbl = value_of(theta_vec[n]);
94 
95  if (include_summand<propto,T_prob>::value) {
96  if (n_int == 1)
97  logp += log(theta_dbl);
98  else
99  logp += log1m(theta_dbl);
100  }
101 
102  // gradient
103  if (include_summand<propto,T_prob>::value) {
104  if (n_int == 1)
105  operands_and_partials.d_x1[n] += 1.0 / theta_dbl;
106  else
107  operands_and_partials.d_x1[n] += 1.0 / (theta_dbl - 1);
108  }
109  }
110  }
111  return operands_and_partials.to_var(logp);
112  }
113 
114  template <bool propto,
115  typename T_y,
116  typename T_prob>
117  inline
118  typename return_type<T_prob>::type
119  bernoulli_log(const T_y& n,
120  const T_prob& theta) {
121  return bernoulli_log<propto>(n,theta,stan::math::default_policy());
122  }
123 
124  template <typename T_y,
125  typename T_prob,
126  class Policy>
127  inline
128  typename return_type<T_prob>::type
129  bernoulli_log(const T_y& n,
130  const T_prob& theta,
131  const Policy&) {
132  return bernoulli_log<false>(n,theta,Policy());
133  }
134 
135  template <typename T_y, typename T_prob>
136  inline
137  typename return_type<T_prob>::type
138  bernoulli_log(const T_y& n,
139  const T_prob& theta) {
140  return bernoulli_log<false>(n,theta,stan::math::default_policy());
141  }
142 
143 
144  // Bernoulli(n|inv_logit(theta)) [0 <= n <= 1; -inf <= theta <= inf]
145  // FIXME: documentation
146  template <bool propto,
147  typename T_n,
148  typename T_prob,
149  class Policy>
150  typename return_type<T_prob>::type
151  bernoulli_logit_log(const T_n& n,
152  const T_prob& theta,
153  const Policy&) {
154  static const char* function = "stan::prob::bernoulli_logit_log(%1%)";
155 
156  using stan::is_constant_struct;
157  using stan::math::check_not_nan;
158  using stan::math::check_bounded;
159  using stan::math::value_of;
160  using stan::math::check_consistent_sizes;
161  using stan::prob::include_summand;
162  using stan::math::log1p;
163  using stan::math::inv_logit;
164 
165  // check if any vectors are zero length
166  if (!(stan::length(n)
167  && stan::length(theta)))
168  return 0.0;
169 
170  // set up return value accumulator
171  double logp(0.0);
172 
173  // validate args (here done over var, which should be OK)
174  if (!check_bounded(function, n, 0, 1, "n", &logp, Policy()))
175  return logp;
176  if (!check_not_nan(function, theta, "Logit transformed probability parameter",
177  &logp, Policy()))
178  return logp;
179  if (!(check_consistent_sizes(function,
180  n,theta,
181  "Random variable","Probability parameter",
182  &logp, Policy())))
183  return logp;
184 
185  // check if no variables are involved and prop-to
186  if (!include_summand<propto,T_prob>::value)
187  return 0.0;
188 
189  // set up template expressions wrapping scalars into vector views
190  VectorView<const T_n> n_vec(n);
191  VectorView<const T_prob> theta_vec(theta);
192  size_t N = max_size(n, theta);
193  agrad::OperandsAndPartials<T_prob> operands_and_partials(theta);
194 
195  for (size_t n = 0; n < N; n++) {
196  // pull out values of arguments
197  const int n_int = value_of(n_vec[n]);
198  const double theta_dbl = value_of(theta_vec[n]);
199 
200  // reusable subexpression values
201  const int sign = 2*n_int-1;
202  const double ntheta = sign * theta_dbl;
203  const double exp_m_ntheta = exp(-ntheta);
204 
205  if (include_summand<propto,T_prob>::value) {
206  // Handle extreme values gracefully using Taylor approximations.
207  const static double cutoff = 20.0;
208  if (ntheta > cutoff)
209  logp -= exp_m_ntheta;
210  else if (ntheta < -cutoff)
211  logp += ntheta;
212  else
213  logp -= log1p(exp_m_ntheta);
214  }
215 
216  // gradients
217  if (!is_constant_struct<T_prob>::value) {
218  const static double cutoff = 20.0;
219  if (ntheta > cutoff)
220  operands_and_partials.d_x1[n] -= exp_m_ntheta;
221  else if (ntheta < -cutoff)
222  operands_and_partials.d_x1[n] += sign;
223  else
224  operands_and_partials.d_x1[n] += sign * exp_m_ntheta / (exp_m_ntheta + 1);
225  }
226  }
227  return operands_and_partials.to_var(logp);
228  }
229 
230 
231  template <bool propto,
232  typename T_n,
233  typename T_prob>
234  inline
235  typename return_type<T_prob>::type
236  bernoulli_logit_log(const T_n& n,
237  const T_prob& theta) {
238  return bernoulli_logit_log<propto>(n,theta,stan::math::default_policy());
239  }
240 
241 
242  template <typename T_n,
243  typename T_prob,
244  class Policy>
245  inline
246  typename return_type<T_prob>::type
247  bernoulli_logit_log(const T_n& n,
248  const T_prob& theta,
249  const Policy&) {
250  return bernoulli_logit_log<false>(n,theta,Policy());
251  }
252 
253 
254  template <typename T_n,
255  typename T_prob>
256  inline
257  typename return_type<T_prob>::type
258  bernoulli_logit_log(const T_n& n,
259  const T_prob& theta) {
260  return bernoulli_logit_log<false>(n,theta,stan::math::default_policy());
261  }
262 
263  // Bernoulli CDF
264  template <bool propto, typename T_n, typename T_prob, class Policy>
265  typename return_type<T_prob>::type
266  bernoulli_cdf(const T_n& n, const T_prob& theta, const Policy&) {
267 
268  static const char* function = "stan::prob::bernoulli_cdf(%1%)";
269 
270  using stan::math::check_finite;
271  using stan::math::check_bounded;
272  using stan::math::check_consistent_sizes;
273  using stan::prob::include_summand;
274 
275  // Ensure non-zero argument lenghts
276  if (!(stan::length(n) && stan::length(theta)))
277  return 0.0;
278 
279  double P(1.0);
280 
281  // Validate arguments
282  if (!check_bounded(function, n, 0, 1, "n", &P, Policy()))
283  return P;
284 
285  if (!check_finite(function, theta, "Probability parameter", &P, Policy()))
286  return P;
287 
288  if (!check_bounded(function, theta, 0.0, 1.0,
289  "Probability parameter", &P, Policy()))
290  return P;
291 
292  if (!(check_consistent_sizes(function,
293  n, theta,
294  "Random variable","Probability parameter",
295  &P, Policy())))
296  return P;
297 
298  // Return if everything is constant and only proportionality is required
299  if (!include_summand<propto,T_prob>::value)
300  return 0.0;
301 
302  // set up template expressions wrapping scalars into vector views
303  VectorView<const T_n> n_vec(n);
304  VectorView<const T_prob> theta_vec(theta);
305  size_t size = max_size(n, theta);
306 
307  // Compute vectorized CDF and gradient
308  using stan::math::value_of;
309  using boost::math::ibetac;
310  using boost::math::ibeta_derivative;
311 
312  agrad::OperandsAndPartials<T_prob> operands_and_partials(theta);
313 
314  for (size_t i = 0; i < size; i++) {
315 
316  const double n_dbl = value_of(n_vec[i]);
317  const double theta_dbl = value_of(theta_vec[i]);
318 
319  const double Pi = ibetac(n_dbl + 1, 1 - n_dbl, theta_dbl);
320 
321  P *= Pi;
322 
323  if (!is_constant_struct<T_prob>::value)
324  operands_and_partials.d_x1[i] += - ibeta_derivative(n_dbl + 1, 1 - n_dbl, theta_dbl) / Pi;
325 
326  }
327 
328  for (size_t i = 0; i < size; i++) {
329 
330  if (!is_constant_struct<T_prob>::value)
331  operands_and_partials.d_x1[i] *= P;
332 
333  }
334 
335  return P;
336 
337  }
338 
339  template <bool propto, typename T_y, typename T_prob>
340  inline typename return_type<T_prob>::type
341  bernoulli_cdf(const T_y& n, const T_prob& theta) {
342  return bernoulli_log<propto>(n, theta, stan::math::default_policy());
343  }
344 
345  template <typename T_y, typename T_prob, class Policy>
346  inline typename return_type<T_prob>::type
347  bernoulli_cdf(const T_y& n, const T_prob& theta, const Policy&) {
348  return bernoulli_cdf<false>(n, theta, Policy());
349  }
350 
351  template <typename T_y, typename T_prob>
352  inline typename return_type<T_prob>::type
353  bernoulli_cdf(const T_y& n, const T_prob& theta) {
354  return bernoulli_cdf<false>(n, theta, stan::math::default_policy());
355  }
356 
357 
358  }
359 }
360 #endif
     [ Stan Home Page ] © 2011–2013, Stan Development Team.