Stan  1.3
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
skew_normal.hpp
Go to the documentation of this file.
1 #ifndef __STAN__PROB__DISTRIBUTIONS__UNIVARIATE__CONTINUOUS__SKEW__NORMAL__HPP__
2 #define __STAN__PROB__DISTRIBUTIONS__UNIVARIATE__CONTINUOUS__SKEW__NORMAL__HPP__
3 
4 #include <boost/random/variate_generator.hpp>
5 #include <boost/math/distributions.hpp>
6 #include <stan/prob/distributions/univariate/continuous/uniform.hpp>
7 
8 #include <stan/agrad.hpp>
9 #include <stan/math/error_handling.hpp>
10 #include <stan/math/functions/owens_t.hpp>
11 #include <stan/agrad/rev/owens_t.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_shape,
23  class Policy>
24  typename return_type<T_y,T_loc,T_scale,T_shape>::type
25  skew_normal_log(const T_y& y, const T_loc& mu, const T_scale& sigma,
26  const T_shape& alpha, const Policy& /*policy*/) {
27  static const char* function = "stan::prob::skew_normal_log(%1%)";
28 
29  using std::log;
30  using stan::is_constant_struct;
31  using stan::math::check_positive;
32  using stan::math::check_finite;
33  using stan::math::check_not_nan;
34  using stan::math::check_consistent_sizes;
35  using stan::math::value_of;
36  using stan::prob::include_summand;
37 
38  // check if any vectors are zero length
39  if (!(stan::length(y)
40  && stan::length(mu)
41  && stan::length(sigma)
42  && stan::length(alpha)))
43  return 0.0;
44 
45  // set up return value accumulator
46  double logp(0.0);
47 
48  // validate args (here done over var, which should be OK)
49  if (!check_not_nan(function, y, "Random variable", &logp, Policy()))
50  return logp;
51  if (!check_finite(function, mu, "Location parameter",
52  &logp, Policy()))
53  return logp;
54  if (!check_finite(function, alpha, "Shape parameter",
55  &logp, Policy()))
56  return logp;
57  if (!check_positive(function, sigma, "Scale parameter",
58  &logp, Policy()))
59  return logp;
60  if (!(check_consistent_sizes(function,
61  y,mu,sigma,alpha,
62  "Random variable","Location parameter","Scale parameter", "Shape paramter",
63  &logp, Policy())))
64  return logp;
65 
66  // check if no variables are involved and prop-to
67  if (!include_summand<propto,T_y,T_loc,T_scale,T_shape>::value)
68  return 0.0;
69 
70  // set up template expressions wrapping scalars into vector views
71  agrad::OperandsAndPartials<T_y, T_loc, T_scale, T_shape> operands_and_partials(y, mu, sigma, alpha);
72 
73  VectorView<const T_y> y_vec(y);
74  VectorView<const T_loc> mu_vec(mu);
75  VectorView<const T_scale> sigma_vec(sigma);
76  VectorView<const T_shape> alpha_vec(alpha);
77  size_t N = max_size(y, mu, sigma, alpha);
78 
79  DoubleVectorView<true,is_vector<T_scale>::value> inv_sigma(length(sigma));
80  DoubleVectorView<include_summand<propto,T_scale>::value,is_vector<T_scale>::value> log_sigma(length(sigma));
81  for (size_t i = 0; i < length(sigma); i++) {
82  inv_sigma[i] = 1.0 / value_of(sigma_vec[i]);
83  if (include_summand<propto,T_scale>::value)
84  log_sigma[i] = log(value_of(sigma_vec[i]));
85  }
86 
87  for (size_t n = 0; n < N; n++) {
88  // pull out values of arguments
89  const double y_dbl = value_of(y_vec[n]);
90  const double mu_dbl = value_of(mu_vec[n]);
91  const double sigma_dbl = value_of(sigma_vec[n]);
92  const double alpha_dbl = value_of(alpha_vec[n]);
93 
94  // reusable subexpression values
95  const double y_minus_mu_over_sigma
96  = (y_dbl - mu_dbl) * inv_sigma[n];
97  const double pi_dbl = boost::math::constants::pi<double>();
98 
99  // log probability
100  if (include_summand<propto>::value)
101  logp -= 0.5 * log(2.0 * pi_dbl);
102  if (include_summand<propto, T_scale>::value)
103  logp -= log(sigma_dbl);
104  if (include_summand<propto,T_y, T_loc, T_scale>::value)
105  logp -= y_minus_mu_over_sigma * y_minus_mu_over_sigma / 2.0;
106  if (include_summand<propto,T_y,T_loc,T_scale,T_shape>::value)
107  logp += log(boost::math::erfc(-alpha_dbl * y_minus_mu_over_sigma / std::sqrt(2.0)));
108 
109  // gradients
110  double deriv_logerf = 2.0 / std::sqrt(pi_dbl) * exp(-alpha_dbl * y_minus_mu_over_sigma / std::sqrt(2.0) * alpha_dbl * y_minus_mu_over_sigma / std::sqrt(2.0)) / (1 + boost::math::erf(alpha_dbl * y_minus_mu_over_sigma / std::sqrt(2.0)));
111  if (!is_constant_struct<T_y>::value)
112  operands_and_partials.d_x1[n] += -y_minus_mu_over_sigma / sigma_dbl + deriv_logerf * alpha_dbl / (sigma_dbl * std::sqrt(2.0)) ;
113  if (!is_constant_struct<T_loc>::value)
114  operands_and_partials.d_x2[n] += y_minus_mu_over_sigma / sigma_dbl + deriv_logerf * -alpha_dbl / (sigma_dbl * std::sqrt(2.0));
115  if (!is_constant_struct<T_scale>::value)
116  operands_and_partials.d_x3[n] += -1.0 / sigma_dbl + y_minus_mu_over_sigma * y_minus_mu_over_sigma / sigma_dbl - deriv_logerf * y_minus_mu_over_sigma * alpha_dbl / (sigma_dbl * std::sqrt(2.0));
117  if (!is_constant_struct<T_shape>::value)
118  operands_and_partials.d_x4[n] += deriv_logerf * y_minus_mu_over_sigma / std::sqrt(2.0);
119  }
120  return operands_and_partials.to_var(logp);
121  }
122 
123 
124  template <bool propto,
125  typename T_y, typename T_loc, typename T_scale, typename T_shape>
126  inline
127  typename return_type<T_y,T_loc,T_scale,T_shape>::type
128  skew_normal_log(const T_y& y, const T_loc& mu, const T_scale& sigma, const T_shape& alpha) {
129  return skew_normal_log<propto>(y,mu,sigma,alpha,stan::math::default_policy());
130  }
131 
132  template <typename T_y, typename T_loc, typename T_scale, typename T_shape,
133  class Policy>
134  inline
135  typename return_type<T_y,T_loc,T_scale,T_shape>::type
136  skew_normal_log(const T_y& y, const T_loc& mu, const T_scale& sigma, const T_shape& alpha, const Policy&) {
137  return skew_normal_log<false>(y,mu,sigma,alpha,Policy());
138  }
139 
140  template <typename T_y, typename T_loc, typename T_scale, typename T_shape>
141  inline
142  typename return_type<T_y,T_loc,T_scale, T_shape>::type
143  skew_normal_log(const T_y& y, const T_loc& mu, const T_scale& sigma, const T_shape& alpha) {
144  return skew_normal_log<false>(y,mu,sigma,alpha,stan::math::default_policy());
145  }
146 
147  template <typename T_y, typename T_loc, typename T_scale, typename T_shape,
148  class Policy>
149  typename return_type<T_y,T_loc,T_scale,T_shape>::type
150  skew_normal_cdf(const T_y& y, const T_loc& mu, const T_scale& sigma, const T_shape& alpha,
151  const Policy&) {
152  static const char* function = "stan::prob::skew_normal_cdf(%1%)";
153 
154  using stan::math::check_positive;
155  using stan::math::check_finite;
156  using stan::math::check_not_nan;
157  using stan::math::check_consistent_sizes;
158  using stan::agrad::owens_t;
159  using stan::math::owens_t;
160 
161 
162  typename return_type<T_y, T_loc, T_scale, T_shape>::type cdf(1);
163  // check if any vectors are zero length
164  if (!(stan::length(y)
165  && stan::length(mu)
166  && stan::length(sigma)
167  && stan::length(alpha)))
168  return cdf;
169 
170  if (!check_not_nan(function, y, "Random variable", &cdf, Policy()))
171  return cdf;
172  if (!check_finite(function, mu, "Location parameter", &cdf, Policy()))
173  return cdf;
174  if (!check_not_nan(function, sigma, "Scale parameter",
175  &cdf, Policy()))
176  return cdf;
177  if (!check_positive(function, sigma, "Scale parameter",
178  &cdf, Policy()))
179  return cdf;
180  if (!check_finite(function, alpha, "Shape parameter", &cdf, Policy()))
181  return cdf;
182  if (!check_not_nan(function, alpha, "Shape parameter",
183  &cdf, Policy()))
184  return cdf;
185  if (!(check_consistent_sizes(function,
186  y,mu,sigma,alpha,
187  "Random variable","Location parameter","Scale parameter","Shape paramter",
188  &cdf, Policy())))
189  return cdf;
190 
191  VectorView<const T_y> y_vec(y);
192  VectorView<const T_loc> mu_vec(mu);
193  VectorView<const T_scale> sigma_vec(sigma);
194  VectorView<const T_shape> alpha_vec(alpha);
195  size_t N = max_size(y, mu, sigma, alpha);
196 
197  for (size_t n = 0; n < N; n++) {
198  cdf *= 0.5 * erfc(-(y_vec[n] - mu_vec[n]) / (std::sqrt(2) * sigma_vec[n])) - 2 * owens_t((y_vec[n] - mu_vec[n]) / sigma_vec[n], alpha_vec[n]);
199  }
200  return cdf;
201  }
202 
203  template <typename T_y, typename T_loc, typename T_scale, typename T_shape>
204  inline
205  typename return_type<T_y, T_loc, T_scale, T_shape>::type
206  skew_normal_cdf(const T_y& y, const T_loc& mu, const T_scale& sigma, const T_shape& alpha) {
207  return skew_normal_cdf(y,mu,sigma,alpha,stan::math::default_policy());
208  }
209 
210  template <class RNG>
211  inline double
212  skew_normal_rng(const double mu,
213  const double sigma,
214  const double alpha,
215  RNG& rng) {
216  boost::math::skew_normal_distribution<>dist (mu, sigma, alpha);
217  return quantile(dist, stan::prob::uniform_rng(0.0,1.0,rng));
218  }
219  }
220 }
221 #endif
222 
     [ Stan Home Page ] © 2011–2013, Stan Development Team.