Stan  1.0
probability, sampling & optimization
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
error_handling.hpp
special_functions.hpp
stan::agrad::log1m
var log1m(const stan::agrad::var &a)
The log (1 - x) function for variables.
Definition: special_functions.hpp:800
stan::agrad::log1p
var log1p(const stan::agrad::var &a)
The log (1 + x) function for variables (C99).
Definition: special_functions.hpp:786
stan::agrad::value_of
double value_of(const agrad::var &v)
Return the value of the specified variable.
Definition: special_functions.hpp:1582
stan::agrad::log
var log(const var &a)
Return the natural log of the specified variable (cmath).
Definition: agrad.hpp:1730
stan::agrad::sum
var sum(const Eigen::Matrix< var, R, C > &m)
Returns the sum of the coefficients of the specified matrix, column vector or row vector.
Definition: matrix.hpp:837
stan::agrad::exp
var exp(const var &a)
Return the exponentiation of the specified variable (cmath).
Definition: agrad.hpp:1716
stan::math::log1m
boost::math::tools::promote_args< T >::type log1m(T x)
Return the natural logarithm of one minus the specified value.
Definition: special_functions.hpp:442
stan::math::log1p
boost::math::tools::promote_args< T >::type log1p(T x)
Return the natural logarithm of one plus the specified value.
Definition: special_functions.hpp:418
stan::math::inv_logit
boost::math::tools::promote_args< T >::type inv_logit(T a)
Returns the inverse logit function applied to the argument.
Definition: special_functions.hpp:205
stan::math::value_of
double value_of(T x)
Return the value of the specified scalar argument converted to a double value.
Definition: special_functions.hpp:842
stan::math::check_not_nan
bool check_not_nan(const char *function, const T_y &y, const char *name, T_result *result, const Policy &)
Checks if the variable y is nan.
Definition: error_handling.hpp:144
stan::math::check_consistent_sizes
bool check_consistent_sizes(const char *function, const T1 &x1, const T2 &x2, const char *name1, const char *name2, T_result *result, const Policy &)
Definition: error_handling.hpp:782
stan::math::check_bounded
bool check_bounded(const char *function, const T_y &y, const T_low &low, const T_high &high, const char *name, T_result *result, const Policy &)
Definition: error_handling.hpp:604
stan::math::check_finite
bool check_finite(const char *function, const T_y &y, const char *name, T_result *result, const Policy &)
Checks if the variable y is finite.
Definition: error_handling.hpp:224
stan::math::default_policy
boost::math::policies::policy default_policy
Default error-handling policy from Boost.
Definition: error_handling.hpp:30
stan::prob::bernoulli_log
return_type< T_prob >::type bernoulli_log(const T_n &n, const T_prob &theta, const Policy &)
Definition: bernoulli.hpp:21
stan::prob::bernoulli_logit_log
return_type< T_prob >::type bernoulli_logit_log(const T_n &n, const T_prob &theta, const Policy &)
Definition: bernoulli.hpp:142
stan
Probability, optimization and sampling library.
Definition: agrad.cpp:6
stan::length
size_t length(const T &x)
Definition: traits.hpp:111
stan::max_size
size_t max_size(const T1 &x1, const T2 &x2)
Definition: traits.hpp:148
stan::agrad::OperandsAndPartials
A variable implementation that stores operands and derivatives with respect to the variable.
Definition: partials_vari.hpp:76
stan::agrad::OperandsAndPartials::d_x1
VectorView< double *, is_vector< T1 >::value > d_x1
Definition: partials_vari.hpp:82
stan::agrad::OperandsAndPartials::to_var
T_return_type to_var(double logp)
Definition: partials_vari.hpp:138
stan::is_constant_struct
Metaprogram to determine if a type has a base scalar type that can be assigned to type double.
Definition: traits.hpp:39
stan::prob::include_summand
Template metaprogram to calculate whether a summand needs to be included in a proportional (log) prob...
Definition: traits.hpp:33
stan::return_type::type
boost::math::tools::promote_args< typename scalar_type< T1 >::type, typename scalar_type< T2 >::type, typename scalar_type< T3 >::type, typename scalar_type< T4 >::type, typename scalar_type< T5 >::type, typename scalar_type< T6 >::type >::type type
Definition: traits.hpp:368
     [ Stan Home Page ] © 2011–2012, Stan Development Team.