Stan  1.3
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
double_exponential.hpp
Go to the documentation of this file.
1 #ifndef __STAN__PROB__DISTRIBUTIONS__DOUBLE_EXPONENTIAL_HPP__
2 #define __STAN__PROB__DISTRIBUTIONS__DOUBLE_EXPONENTIAL_HPP__
3 
4 #include <boost/random/uniform_01.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/value_of.hpp>
10 #include <stan/meta/traits.hpp>
11 #include <stan/prob/constants.hpp>
12 #include <stan/prob/traits.hpp>
13 
14 namespace stan {
15 
16  namespace prob {
17  template<typename T>
18  inline int sign(const T& z) {
19  return (z == 0) ? 0 : z < 0 ? -1 : 1;
20  }
21 
22  // DoubleExponential(y|mu,sigma) [sigma > 0]
23  // FIXME: add documentation
24  template <bool propto,
25  typename T_y, typename T_loc, typename T_scale,
26  class Policy>
27  typename return_type<T_y,T_loc,T_scale>::type
28  double_exponential_log(const T_y& y, const T_loc& mu, const T_scale& sigma,
29  const Policy&) {
30  static const char* function
31  = "stan::prob::double_exponential_log(%1%)";
32 
33  using stan::is_constant_struct;
34  using stan::math::check_finite;
35  using stan::math::check_positive;
36  using stan::math::check_consistent_sizes;
37  using stan::math::value_of;
38  using stan::prob::include_summand;
39  using std::log;
40  using std::fabs;
41 
42  // check if any vectors are zero length
43  if (!(stan::length(y)
44  && stan::length(mu)
45  && stan::length(sigma)))
46  return 0.0;
47 
48  // set up return value accumulator
49  double logp(0.0);
50  if(!check_finite(function, y, "Random variable", &logp, Policy()))
51  return logp;
52  if(!check_finite(function, mu, "Location parameter",
53  &logp, Policy()))
54  return logp;
55  if(!check_finite(function, sigma, "Scale parameter",
56  &logp, Policy()))
57  return logp;
58  if(!check_positive(function, sigma, "Scale parameter",
59  &logp, Policy()))
60  return logp;
61  if (!(check_consistent_sizes(function,
62  y,mu,sigma,
63  "Random variable","Location parameter","Shape parameter",
64  &logp, Policy())))
65  return logp;
66 
67  // check if no variables are involved and prop-to
68  if (!include_summand<propto,T_y,T_loc,T_scale>::value)
69  return 0.0;
70 
71  // set up template expressions wrapping scalars into vector views
72  VectorView<const T_y> y_vec(y);
73  VectorView<const T_loc> mu_vec(mu);
74  VectorView<const T_scale> sigma_vec(sigma);
75  size_t N = max_size(y, mu, sigma);
76  agrad::OperandsAndPartials<T_y,T_loc,T_scale> operands_and_partials(y, mu, sigma);
77 
78  DoubleVectorView<include_summand<propto,T_y,T_loc,T_scale>::value,is_vector<T_scale>::value>
79  inv_sigma(length(sigma));
80  DoubleVectorView<!is_constant_struct<T_scale>::value,is_vector<T_scale>::value>
81  inv_sigma_squared(length(sigma));
82  DoubleVectorView<include_summand<propto,T_scale>::value,is_vector<T_scale>::value>
83  log_sigma(length(sigma));
84  for (size_t i = 0; i < length(sigma); i++) {
85  const double sigma_dbl = value_of(sigma_vec[i]);
86  if (include_summand<propto,T_y,T_loc,T_scale>::value)
87  inv_sigma[i] = 1.0 / sigma_dbl;
88  if (include_summand<propto,T_scale>::value)
89  log_sigma[i] = log(value_of(sigma_vec[i]));
90  if (!is_constant_struct<T_scale>::value)
91  inv_sigma_squared[i] = inv_sigma[i] * inv_sigma[i];
92  }
93 
94 
95  for (size_t n = 0; n < N; n++) {
96  const double y_dbl = value_of(y_vec[n]);
97  const double mu_dbl = value_of(mu_vec[n]);
98 
99  // reusable subexpressions values
100  const double y_m_mu = y_dbl - mu_dbl;
101  const double fabs_y_m_mu = fabs(y_m_mu);
102 
103  // log probability
104  if (include_summand<propto>::value)
105  logp += NEG_LOG_TWO;
106  if (include_summand<propto,T_scale>::value)
107  logp -= log_sigma[n];
108  if (include_summand<propto,T_y,T_loc,T_scale>::value)
109  logp -= fabs_y_m_mu * inv_sigma[n];
110 
111  // gradients
112  double sign_y_m_mu_times_inv_sigma(0);
113  if (!is_constant_struct<T_y>::value || !is_constant_struct<T_loc>::value)
114  sign_y_m_mu_times_inv_sigma = sign(y_m_mu) * inv_sigma[n];
115  if (!is_constant_struct<T_y>::value) {
116  operands_and_partials.d_x1[n] -= sign_y_m_mu_times_inv_sigma;
117  }
118  if (!is_constant_struct<T_loc>::value) {
119  operands_and_partials.d_x2[n] += sign_y_m_mu_times_inv_sigma;
120  }
121  if (!is_constant_struct<T_scale>::value)
122  operands_and_partials.d_x3[n] += -inv_sigma[n] + fabs_y_m_mu * inv_sigma_squared[n];
123  }
124  return operands_and_partials.to_var(logp);
125  }
126 
127 
128  template <bool propto,
129  typename T_y, typename T_loc, typename T_scale>
130  typename return_type<T_y,T_loc,T_scale>::type
131  double_exponential_log(const T_y& y, const T_loc& mu,
132  const T_scale& sigma) {
133  return double_exponential_log<propto>(y,mu,sigma,
134  stan::math::default_policy());
135  }
136 
137 
138  template <typename T_y, typename T_loc, typename T_scale,
139  class Policy>
140  typename return_type<T_y,T_loc,T_scale>::type
141  double_exponential_log(const T_y& y, const T_loc& mu, const T_scale& sigma,
142  const Policy&) {
143  return double_exponential_log<false>(y,mu,sigma,Policy());
144  }
145 
146  template <typename T_y, typename T_loc, typename T_scale>
147  typename return_type<T_y,T_loc,T_scale>::type
148  double_exponential_log(const T_y& y, const T_loc& mu,
149  const T_scale& sigma) {
150  return double_exponential_log<false>(y,mu,sigma,
151  stan::math::default_policy());
152  }
153 
168  template <typename T_y, typename T_loc, typename T_scale,
169  class Policy>
170  typename return_type<T_y,T_loc,T_scale>::type
171  double_exponential_cdf(const T_y& y, const T_loc& mu, const T_scale& sigma,
172  const Policy&) {
173  static const char* function
174  = "stan::prob::double_exponential_cdf(%1%)";
175 
176  using stan::math::check_finite;
177  using stan::math::check_positive;
178  using boost::math::tools::promote_args;
179 
180  typename promote_args<T_y,T_loc,T_scale>::type lp(0.0);
181  if(!check_finite(function, y, "Random variable", &lp, Policy()))
182  return lp;
183  if(!check_finite(function, mu, "Location parameter",
184  &lp, Policy()))
185  return lp;
186  if(!check_finite(function, sigma, "Scale parameter",
187  &lp, Policy()))
188  return lp;
189  if(!check_positive(function, sigma, "Scale parameter",
190  &lp, Policy()))
191  return lp;
192 
193  if (y < mu)
194  return exp((y-mu)/sigma)/2;
195  else
196  return 1 - exp((mu-y)/sigma)/2;
197  }
198 
199  template <typename T_y, typename T_loc, typename T_scale>
200  typename boost::math::tools::promote_args<T_y,T_loc,T_scale>::type
201  double_exponential_cdf(const T_y& y, const T_loc& mu, const T_scale& sigma) {
202  return double_exponential_cdf(y,mu,sigma,stan::math::default_policy());
203  }
204 
205  template <class RNG>
206  inline double
207  double_exponential_rng(const double mu,
208  const double sigma,
209  RNG& rng) {
210  using boost::variate_generator;
211  using boost::random::uniform_01;
212  using std::log;
213  using std::abs;
214  variate_generator<RNG&, uniform_01<> >
215  rng_unit_01(rng, uniform_01<>());
216  double a = 0;
217  double laplaceRN = rng_unit_01();
218  if(0.5 - laplaceRN > 0)
219  a = 1.0;
220  else if(0.5 - laplaceRN < 0)
221  a = -1.0;
222  return mu - sigma * a * log(1 - 2 * abs(0.5 - laplaceRN));
223  }
224  }
225 }
226 #endif
     [ Stan Home Page ] © 2011–2013, Stan Development Team.