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