Stan  1.3
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
exp_mod_normal.hpp
Go to the documentation of this file.
1 #ifndef __STAN__PROB__DISTRIBUTIONS__UNIVARIATE__CONTINUOUS__EXP__MOD__NORMAL__HPP__
2 #define __STAN__PROB__DISTRIBUTIONS__UNIVARIATE__CONTINUOUS__EXP__MOD__NORMAL__HPP__
3 
4 #include <boost/random/normal_distribution.hpp>
5 #include <boost/math/special_functions/fpclassify.hpp>
6 #include <boost/random/variate_generator.hpp>
7 #include <stan/prob/distributions/univariate/continuous/normal.hpp>
8 #include <stan/prob/distributions/univariate/continuous/exponential.hpp>
9 
10 #include <stan/agrad.hpp>
11 #include <stan/math/error_handling.hpp>
12 #include <stan/meta/traits.hpp>
13 #include <stan/prob/constants.hpp>
14 #include <stan/prob/traits.hpp>
15 #include <stan/math/functions/value_of.hpp>
16 
17 namespace stan {
18 
19  namespace prob {
20 
21  template <bool propto,
22  typename T_y, typename T_loc, typename T_scale,typename T_inv_scale, class Policy>
23  typename return_type<T_y,T_loc,T_scale, T_inv_scale>::type
24  exp_mod_normal_log(const T_y& y, const T_loc& mu, const T_scale& sigma,
25  const T_inv_scale& lambda, const Policy& /*policy*/) {
26  static const char* function = "stan::prob::exp_mod_normal_log(%1%)";
27 
28  using stan::is_constant_struct;
29  using stan::math::check_positive;
30  using stan::math::check_finite;
31  using stan::math::check_not_nan;
32  using stan::math::check_consistent_sizes;
33  using stan::math::value_of;
34  using stan::prob::include_summand;
35 
36  // check if any vectors are zero length
37  if (!(stan::length(y)
38  && stan::length(mu)
39  && stan::length(sigma)
40  && stan::length(lambda)))
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_finite(function, lambda, "Inv_scale parameter",
53  &logp, Policy()))
54  return logp;
55  if (!check_positive(function, lambda, "Inv_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,lambda,
63  "Random variable","Location parameter","Scale parameter", "Inv_scale paramter",
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,T_inv_scale>::value)
69  return 0.0;
70 
71  // set up template expressions wrapping scalars into vector views
72  agrad::OperandsAndPartials<T_y, T_loc, T_scale, T_inv_scale> operands_and_partials(y, mu, sigma,lambda);
73 
74  VectorView<const T_y> y_vec(y);
75  VectorView<const T_loc> mu_vec(mu);
76  VectorView<const T_scale> sigma_vec(sigma);
77  VectorView<const T_inv_scale> lambda_vec(lambda);
78  size_t N = max_size(y, mu, sigma, lambda);
79 
80  for (size_t n = 0; n < N; n++) {
81  //pull out values of arguments
82  const double y_dbl = value_of(y_vec[n]);
83  const double mu_dbl = value_of(mu_vec[n]);
84  const double sigma_dbl = value_of(sigma_vec[n]);
85  const double lambda_dbl = value_of(lambda_vec[n]);
86 
87  const double pi_dbl = boost::math::constants::pi<double>();
88 
89  // log probability
90  if (include_summand<propto>::value)
91  logp -= log(2.0);
92  if (include_summand<propto, T_inv_scale>::value)
93  logp += log(lambda_dbl);
94  if (include_summand<propto,T_y,T_loc,T_scale,T_inv_scale>::value)
95  logp += lambda_dbl * (mu_dbl + 0.5 * lambda_dbl * sigma_dbl * sigma_dbl - y_dbl) + log(boost::math::erfc((mu_dbl + lambda_dbl * sigma_dbl * sigma_dbl - y_dbl) / (std::sqrt(2.0) * sigma_dbl)));
96 
97  // gradients
98  const double deriv_logerfc = -2.0 / std::sqrt(pi_dbl) * exp(-(mu_dbl + lambda_dbl * sigma_dbl * sigma_dbl - y_dbl) / (std::sqrt(2.0) * sigma_dbl) * (mu_dbl + lambda_dbl * sigma_dbl * sigma_dbl - y_dbl) / (sigma_dbl * std::sqrt(2.0))) / boost::math::erfc((mu_dbl + lambda_dbl * sigma_dbl * sigma_dbl - y_dbl) / (sigma_dbl * std::sqrt(2.0)));
99 
100  if (!is_constant_struct<T_y>::value)
101  operands_and_partials.d_x1[n] += -lambda_dbl + deriv_logerfc * -1.0 / (sigma_dbl * std::sqrt(2.0));
102  if (!is_constant_struct<T_loc>::value)
103  operands_and_partials.d_x2[n] += lambda_dbl + deriv_logerfc / (sigma_dbl * std::sqrt(2.0));
104  if (!is_constant_struct<T_scale>::value)
105  operands_and_partials.d_x3[n] += sigma_dbl * lambda_dbl * lambda_dbl + deriv_logerfc * (-mu_dbl / (sigma_dbl * sigma_dbl * std::sqrt(2.0)) + lambda_dbl / std::sqrt(2.0) + y_dbl / (sigma_dbl * sigma_dbl * std::sqrt(2.0)));
106  if (!is_constant_struct<T_inv_scale>::value)
107  operands_and_partials.d_x4[n] += 1 / lambda_dbl + lambda_dbl * sigma_dbl * sigma_dbl + mu_dbl - y_dbl + deriv_logerfc * sigma_dbl / std::sqrt(2.0);
108  }
109  return operands_and_partials.to_var(logp);
110  }
111 
112  template <bool propto,
113  typename T_y, typename T_loc, typename T_scale, typename T_inv_scale>
114  inline
115  typename return_type<T_y,T_loc,T_scale, T_inv_scale>::type
116  exp_mod_normal_log(const T_y& y, const T_loc& mu, const T_scale& sigma, const T_inv_scale& lambda) {
117  return exp_mod_normal_log<propto>(y,mu,sigma,lambda,stan::math::default_policy());
118  }
119 
120  template <typename T_y, typename T_loc, typename T_scale, typename T_inv_scale,
121  class Policy>
122  inline
123  typename return_type<T_y,T_loc,T_scale, T_inv_scale>::type
124  exp_mod_normal_log(const T_y& y, const T_loc& mu, const T_scale& sigma, const T_inv_scale& lambda, const Policy&) {
125  return exp_mod_normal_log<false>(y,mu,sigma,lambda,Policy());
126  }
127 
128  template <typename T_y, typename T_loc, typename T_scale, typename T_inv_scale>
129  inline
130  typename return_type<T_y,T_loc,T_scale, T_inv_scale>::type
131  exp_mod_normal_log(const T_y& y, const T_loc& mu, const T_scale& sigma, const T_inv_scale& lambda) {
132  return exp_mod_normal_log<false>(y,mu,sigma,lambda,stan::math::default_policy());
133  }
134 
135  template <typename T_y, typename T_loc, typename T_scale, typename T_inv_scale,
136  class Policy>
137  typename return_type<T_y,T_loc,T_scale,T_inv_scale>::type
138  exp_mod_normal_cdf(const T_y& y, const T_loc& mu, const T_scale& sigma, const T_inv_scale& lambda,
139  const Policy&) {
140  static const char* function = "stan::prob::exp_mod_normal_cdf(%1%)";
141 
142  using stan::math::check_positive;
143  using stan::math::check_finite;
144  using stan::math::check_not_nan;
145  using stan::math::check_consistent_sizes;
146 
147  typename return_type<T_y, T_loc, T_scale, T_inv_scale>::type cdf(1);
148 
149  //check if any vectors are zero length
150  if (!(stan::length(y)
151  && stan::length(mu)
152  && stan::length(sigma)
153  && stan::length(lambda)))
154  return cdf;
155 
156  if (!check_not_nan(function, y, "Random variable", &cdf, Policy()))
157  return cdf;
158  if (!check_finite(function, mu, "Location parameter", &cdf, Policy()))
159  return cdf;
160  if (!check_not_nan(function, sigma, "Scale parameter",
161  &cdf, Policy()))
162  return cdf;
163  if (!check_finite(function, sigma, "Scale parameter", &cdf, Policy()))
164  return cdf;
165  if (!check_positive(function, sigma, "Scale parameter",
166  &cdf, Policy()))
167  return cdf;
168  if (!check_finite(function, lambda, "Inv_scale parameter", &cdf, Policy()))
169  return cdf;
170  if (!check_positive(function, lambda, "Inv_scale parameter",
171  &cdf, Policy()))
172  return cdf;
173  if (!(check_consistent_sizes(function,
174  y,mu,sigma,lambda,
175  "Random variable","Location parameter","Scale parameter","Inv_scale paramter",
176  &cdf, Policy())))
177  return cdf;
178 
179  VectorView<const T_y> y_vec(y);
180  VectorView<const T_loc> mu_vec(mu);
181  VectorView<const T_scale> sigma_vec(sigma);
182  VectorView<const T_inv_scale> lambda_vec(lambda);
183  size_t N = max_size(y, mu, sigma, lambda);
184 
185  for (size_t n = 0; n < N; n++) {
186  if(boost::math::isinf(y_vec[n]))
187  {
188  if (y_vec[n] < 0.0)
189  return cdf * 0.0;
190  }
191  cdf *= 0.5 * (1 + erf((y_vec[n] - mu_vec[n]) / (sqrt(2.0) * sigma_vec[n]))) - exp(-lambda_vec[n] * (y_vec[n] - mu_vec[n]) + lambda_vec[n] * sigma_vec[n] * lambda_vec[n] * sigma_vec[n] / 2.0) * (0.5 * (1 + erf((y_vec[n] - mu_vec[n] - sigma_vec[n] * lambda_vec[n] * sigma_vec[n]) / (sqrt(2.0) * sigma_vec[n]))));
192  }
193 
194  return cdf;
195  }
196 
197  template <typename T_y, typename T_loc, typename T_scale, typename T_inv_scale>
198  inline
199  typename return_type<T_y, T_loc, T_scale, T_inv_scale>::type
200  exp_mod_normal_cdf(const T_y& y, const T_loc& mu, const T_scale& sigma, const T_inv_scale& lambda) {
201  return exp_mod_normal_cdf(y,mu,sigma,lambda,stan::math::default_policy());
202  }
203 
204  template <class RNG>
205  inline double
206  exp_mod_normal_rng(const double mu,
207  const double sigma,
208  const double lambda,
209  RNG& rng) {
210  return stan::prob::normal_rng(mu, sigma,rng) + stan::prob::exponential_rng(lambda, rng);
211  }
212  }
213 }
214 #endif
215 
216 
217 
     [ Stan Home Page ] © 2011–2013, Stan Development Team.