Stan  1.0
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  const double log_theta = log(theta_dbl);
72  const double log1m_theta = log1m(theta_dbl);
73  if (include_summand<propto,T_prob>::value) {
74  logp += sum * log_theta;
75  logp += (N - sum) * log1m_theta;
76 
77  operands_and_partials.d_x1[0] += sum / theta_dbl;
78  operands_and_partials.d_x1[0] += (N - sum) / (theta_dbl - 1);
79  }
80  } else {
81  for (size_t n = 0; n < N; n++) {
82  // pull out values of arguments
83  const int n_int = value_of(n_vec[n]);
84  const double theta_dbl = value_of(theta_vec[n]);
85 
86  if (include_summand<propto,T_prob>::value) {
87  if (n_int == 1)
88  logp += log(theta_dbl);
89  else
90  logp += log1m(theta_dbl);
91  }
92 
93  // gradient
94  if (include_summand<propto,T_prob>::value) {
95  if (n_int == 1)
96  operands_and_partials.d_x1[n] += 1.0 / theta_dbl;
97  else
98  operands_and_partials.d_x1[n] += 1.0 / (theta_dbl - 1);
99  }
100  }
101  }
102  return operands_and_partials.to_var(logp);
103  }
104 
105  template <bool propto,
106  typename T_y,
107  typename T_prob>
108  inline
109  typename return_type<T_prob>::type
110  bernoulli_log(const T_y& n,
111  const T_prob& theta) {
112  return bernoulli_log<propto>(n,theta,stan::math::default_policy());
113  }
114 
115  template <typename T_y,
116  typename T_prob,
117  class Policy>
118  inline
119  typename return_type<T_prob>::type
120  bernoulli_log(const T_y& n,
121  const T_prob& theta,
122  const Policy&) {
123  return bernoulli_log<false>(n,theta,Policy());
124  }
125 
126  template <typename T_y, typename T_prob>
127  inline
128  typename return_type<T_prob>::type
129  bernoulli_log(const T_y& n,
130  const T_prob& theta) {
131  return bernoulli_log<false>(n,theta,stan::math::default_policy());
132  }
133 
134 
135  // Bernoulli(n|inv_logit(theta)) [0 <= n <= 1; -inf <= theta <= inf]
136  // FIXME: documentation
137  template <bool propto,
138  typename T_n,
139  typename T_prob,
140  class Policy>
141  typename return_type<T_prob>::type
142  bernoulli_logit_log(const T_n& n,
143  const T_prob& theta,
144  const Policy&) {
145  static const char* function = "stan::prob::bernoulli_logit_log(%1%)";
146 
147  using stan::is_constant_struct;
148  using stan::math::check_not_nan;
149  using stan::math::check_bounded;
150  using stan::math::value_of;
151  using stan::math::check_consistent_sizes;
152  using stan::prob::include_summand;
153  using stan::math::log1p;
154  using stan::math::inv_logit;
155 
156  // check if any vectors are zero length
157  if (!(stan::length(n)
158  && stan::length(theta)))
159  return 0.0;
160 
161  // set up return value accumulator
162  double logp(0.0);
163 
164  // validate args (here done over var, which should be OK)
165  if (!check_bounded(function, n, 0, 1, "n", &logp, Policy()))
166  return logp;
167  if (!check_not_nan(function, theta, "Logit transformed probability parameter",
168  &logp, Policy()))
169  return logp;
170  if (!(check_consistent_sizes(function,
171  n,theta,
172  "Random variable","Probability parameter",
173  &logp, Policy())))
174  return logp;
175 
176  // check if no variables are involved and prop-to
177  if (!include_summand<propto,T_prob>::value)
178  return 0.0;
179 
180  // set up template expressions wrapping scalars into vector views
181  VectorView<const T_n> n_vec(n);
182  VectorView<const T_prob> theta_vec(theta);
183  size_t N = max_size(n, theta);
184  agrad::OperandsAndPartials<T_prob> operands_and_partials(theta);
185 
186  for (size_t n = 0; n < N; n++) {
187  // pull out values of arguments
188  const int n_int = value_of(n_vec[n]);
189  const double theta_dbl = value_of(theta_vec[n]);
190 
191  // reusable subexpression values
192  const int sign = 2*n_int-1;
193  const double ntheta = sign * theta_dbl;
194  const double exp_m_ntheta = exp(-ntheta);
195 
196  if (include_summand<propto,T_prob>::value) {
197  // Handle extreme values gracefully using Taylor approximations.
198  const static double cutoff = 20.0;
199  if (ntheta > cutoff)
200  logp -= exp_m_ntheta;
201  else if (ntheta < -cutoff)
202  logp += ntheta;
203  else
204  logp -= log1p(exp_m_ntheta);
205  }
206 
207  // gradients
208  if (!is_constant_struct<T_prob>::value) {
209  const static double cutoff = 20.0;
210  if (ntheta > cutoff)
211  operands_and_partials.d_x1[n] -= exp_m_ntheta;
212  else if (ntheta < -cutoff)
213  operands_and_partials.d_x1[n] += sign;
214  else
215  operands_and_partials.d_x1[n] += sign * exp_m_ntheta / (exp_m_ntheta + 1);
216  }
217  }
218  return operands_and_partials.to_var(logp);
219  }
220 
221 
222  template <bool propto,
223  typename T_n,
224  typename T_prob>
225  inline
226  typename return_type<T_prob>::type
227  bernoulli_logit_log(const T_n& n,
228  const T_prob& theta) {
229  return bernoulli_logit_log<propto>(n,theta,stan::math::default_policy());
230  }
231 
232 
233  template <typename T_n,
234  typename T_prob,
235  class Policy>
236  inline
237  typename return_type<T_prob>::type
238  bernoulli_logit_log(const T_n& n,
239  const T_prob& theta,
240  const Policy&) {
241  return bernoulli_logit_log<false>(n,theta,Policy());
242  }
243 
244 
245  template <typename T_n,
246  typename T_prob>
247  inline
248  typename return_type<T_prob>::type
249  bernoulli_logit_log(const T_n& n,
250  const T_prob& theta) {
251  return bernoulli_logit_log<false>(n,theta,stan::math::default_policy());
252  }
253 
254  }
255 }
256 #endif
     [ Stan Home Page ] © 2011–2012, Stan Development Team.