Stan  1.3
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
gumbel.hpp
Go to the documentation of this file.
1 #ifndef __STAN__PROB__DISTRIBUTIONS__UNIVARIATE__CONTINUOUS__GUMBEL_HPP__
2 #define __STAN__PROB__DISTRIBUTIONS__UNIVARIATE__CONTINUOUS__GUMBEL_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/meta/traits.hpp>
10 #include <stan/prob/constants.hpp>
11 #include <stan/prob/traits.hpp>
12 #include <stan/prob/internal_math.hpp>
13 #include <stan/math/functions/value_of.hpp>
14 
15 namespace stan {
16 
17  namespace prob {
18 
19  template <bool propto,
20  typename T_y, typename T_loc, typename T_scale,
21  class Policy>
22  typename return_type<T_y,T_loc,T_scale>::type
23  gumbel_log(const T_y& y, const T_loc& mu, const T_scale& beta,
24  const Policy& /*policy*/) {
25  static const char* function = "stan::prob::gumbel_log(%1%)";
26 
27  using std::log;
28  using std::exp;
29  using stan::is_constant_struct;
30  using stan::math::check_positive;
31  using stan::math::check_finite;
32  using stan::math::check_not_nan;
33  using stan::math::check_consistent_sizes;
34  using stan::math::value_of;
35  using stan::prob::include_summand;
36 
37  // check if any vectors are zero length
38  if (!(stan::length(y)
39  && stan::length(mu)
40  && stan::length(beta)))
41  return 0.0;
42 
43  // set up return value accumulator
44  double logp(0.0);
45 
46  // validate args (here done over var, which should be OK)
47  if (!check_not_nan(function, y, "Random variable", &logp, Policy()))
48  return logp;
49  if (!check_finite(function, mu, "Location parameter",
50  &logp, Policy()))
51  return logp;
52  if (!check_positive(function, beta, "Scale parameter",
53  &logp, Policy()))
54  return logp;
55  if (!(check_consistent_sizes(function,
56  y,mu,beta,
57  "Random variable","Location parameter","Scale parameter",
58  &logp, Policy())))
59  return logp;
60 
61  // check if no variables are involved and prop-to
62  if (!include_summand<propto,T_y,T_loc,T_scale>::value)
63  return 0.0;
64 
65  // set up template expressions wrapping scalars into vector views
66  agrad::OperandsAndPartials<T_y, T_loc, T_scale> operands_and_partials(y, mu, beta);
67 
68  VectorView<const T_y> y_vec(y);
69  VectorView<const T_loc> mu_vec(mu);
70  VectorView<const T_scale> beta_vec(beta);
71  size_t N = max_size(y, mu, beta);
72 
73  DoubleVectorView<true,is_vector<T_scale>::value> inv_beta(length(beta));
74  DoubleVectorView<include_summand<propto,T_scale>::value,is_vector<T_scale>::value> log_beta(length(beta));
75  for (size_t i = 0; i < length(beta); i++) {
76  inv_beta[i] = 1.0 / value_of(beta_vec[i]);
77  if (include_summand<propto,T_scale>::value)
78  log_beta[i] = log(value_of(beta_vec[i]));
79  }
80 
81  for (size_t n = 0; n < N; n++) {
82  // pull out values of arguments
83  const double y_dbl = value_of(y_vec[n]);
84  const double mu_dbl = value_of(mu_vec[n]);
85 
86  // reusable subexpression values
87  const double y_minus_mu_over_beta
88  = (y_dbl - mu_dbl) * inv_beta[n];
89 
90  // log probability
91  if (include_summand<propto,T_scale>::value)
92  logp -= log_beta[n];
93  if (include_summand<propto,T_y,T_loc,T_scale>::value)
94  logp += -y_minus_mu_over_beta - exp(-y_minus_mu_over_beta);
95 
96  // gradients
97  double scaled_diff = inv_beta[n] * exp(-y_minus_mu_over_beta);
98  if (!is_constant_struct<T_y>::value)
99  operands_and_partials.d_x1[n] -= inv_beta[n] - scaled_diff;
100  if (!is_constant_struct<T_loc>::value)
101  operands_and_partials.d_x2[n] += inv_beta[n] - scaled_diff;
102  if (!is_constant_struct<T_scale>::value)
103  operands_and_partials.d_x3[n]
104  += -inv_beta[n] + y_minus_mu_over_beta * inv_beta[n] - scaled_diff * y_minus_mu_over_beta;
105  }
106  return operands_and_partials.to_var(logp);
107  }
108 
109  template <bool propto,
110  typename T_y, typename T_loc, typename T_scale>
111  inline
112  typename return_type<T_y,T_loc,T_scale>::type
113  gumbel_log(const T_y& y, const T_loc& mu, const T_scale& beta) {
114  return gumbel_log<propto>(y,mu,beta,stan::math::default_policy());
115  }
116 
117  template <typename T_y, typename T_loc, typename T_scale,
118  class Policy>
119  inline
120  typename return_type<T_y,T_loc,T_scale>::type
121  gumbel_log(const T_y& y, const T_loc& mu, const T_scale& beta,
122  const Policy&) {
123  return gumbel_log<false>(y,mu,beta,Policy());
124  }
125 
126  template <typename T_y, typename T_loc, typename T_scale>
127  inline
128  typename return_type<T_y,T_loc,T_scale>::type
129  gumbel_log(const T_y& y, const T_loc& mu, const T_scale& beta) {
130  return gumbel_log<false>(y,mu,beta,stan::math::default_policy());
131  }
132 
133  template <typename T_y, typename T_loc, typename T_scale,
134  class Policy>
135  typename return_type<T_y,T_loc,T_scale>::type
136  gumbel_cdf(const T_y& y, const T_loc& mu, const T_scale& beta,
137  const Policy&) {
138  static const char* function = "stan::prob::gumbel_cdf(%1%)";
139 
140  using stan::math::check_positive;
141  using stan::math::check_finite;
142  using stan::math::check_not_nan;
143  using stan::math::check_consistent_sizes;
144 
145  typename return_type<T_y, T_loc, T_scale>::type cdf(1);
146  // check if any vectors are zero length
147  if (!(stan::length(y)
148  && stan::length(mu)
149  && stan::length(beta)))
150  return cdf;
151 
152  if (!check_not_nan(function, y, "Random variable", &cdf, Policy()))
153  return cdf;
154  if (!check_finite(function, mu, "Location parameter", &cdf, Policy()))
155  return cdf;
156  if (!check_not_nan(function, beta, "Scale parameter",
157  &cdf, Policy()))
158  return cdf;
159  if (!check_positive(function, beta, "Scale parameter",
160  &cdf, Policy()))
161  return cdf;
162  if (!(check_consistent_sizes(function,
163  y,mu,beta,
164  "Random variable","Location parameter","Scale parameter",
165  &cdf, Policy())))
166  return cdf;
167 
168  VectorView<const T_y> y_vec(y);
169  VectorView<const T_loc> mu_vec(mu);
170  VectorView<const T_scale> beta_vec(beta);
171  size_t N = max_size(y, mu, beta);
172 
173  for (size_t n = 0; n < N; n++) {
174  cdf *= exp(-exp(-((y_vec[n]) - (mu_vec[n])) / (beta_vec[n])));
175  }
176 
177  return cdf;
178  }
179 
180  template <typename T_y, typename T_loc, typename T_scale>
181  inline
182  typename return_type<T_y, T_loc, T_scale>::type
183  gumbel_cdf(const T_y& y, const T_loc& mu, const T_scale& beta) {
184  return gumbel_cdf(y,mu,beta,stan::math::default_policy());
185  }
186 
187 
188  template <class RNG>
189  inline double
190  gumbel_rng(const double mu,
191  const double beta,
192  RNG& rng) {
193  using boost::variate_generator;
194  using boost::uniform_01;
195  variate_generator<RNG&, uniform_01<> >
196  uniform01_rng(rng, uniform_01<>());
197  return mu - beta * std::log(-std::log(uniform01_rng()));
198  }
199  }
200 }
201 #endif
202 
     [ Stan Home Page ] © 2011–2013, Stan Development Team.