Stan  1.0
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
normal.hpp
Go to the documentation of this file.
1 #ifndef __STAN__PROB__DISTRIBUTIONS__UNIVARIATE__CONTINUOUS__NORMAL_HPP__
2 #define __STAN__PROB__DISTRIBUTIONS__UNIVARIATE__CONTINUOUS__NORMAL_HPP__
3 
4 #include <boost/random/normal_distribution.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/special_functions.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 
36  template <bool propto,
37  typename T_y, typename T_loc, typename T_scale,
38  class Policy>
39  typename return_type<T_y,T_loc,T_scale>::type
40  normal_log(const T_y& y, const T_loc& mu, const T_scale& sigma,
41  const Policy& /*policy*/) {
42  static const char* function = "stan::prob::normal_log(%1%)";
43 
44  using std::log;
45  using stan::is_constant_struct;
46  using stan::math::check_positive;
47  using stan::math::check_finite;
48  using stan::math::check_not_nan;
49  using stan::math::check_consistent_sizes;
50  using stan::math::value_of;
51  using stan::prob::include_summand;
52 
53  // check if any vectors are zero length
54  if (!(stan::length(y)
55  && stan::length(mu)
56  && stan::length(sigma)))
57  return 0.0;
58 
59  // set up return value accumulator
60  double logp(0.0);
61 
62  // validate args (here done over var, which should be OK)
63  if (!check_not_nan(function, y, "Random variable", &logp, Policy()))
64  return logp;
65  if (!check_finite(function, mu, "Location parameter",
66  &logp, Policy()))
67  return logp;
68  if (!check_positive(function, sigma, "Scale parameter",
69  &logp, Policy()))
70  return logp;
71  if (!(check_consistent_sizes(function,
72  y,mu,sigma,
73  "Random variable","Location parameter","Scale parameter",
74  &logp, Policy())))
75  return logp;
76 
77  // check if no variables are involved and prop-to
78  if (!include_summand<propto,T_y,T_loc,T_scale>::value)
79  return 0.0;
80 
81  // set up template expressions wrapping scalars into vector views
82  agrad::OperandsAndPartials<T_y, T_loc, T_scale> operands_and_partials(y, mu, sigma);
83 
84  VectorView<const T_y> y_vec(y);
85  VectorView<const T_loc> mu_vec(mu);
86  VectorView<const T_scale> sigma_vec(sigma);
87  size_t N = max_size(y, mu, sigma);
88 
89  DoubleVectorView<true,T_scale> inv_sigma(length(sigma));
90  DoubleVectorView<include_summand<propto,T_scale>::value,T_scale> log_sigma(length(sigma));
91  for (size_t i = 0; i < length(sigma); i++) {
92  inv_sigma[i] = 1.0 / value_of(sigma_vec[i]);
93  if (include_summand<propto,T_scale>::value)
94  log_sigma[i] = log(value_of(sigma_vec[i]));
95  }
96 
97  for (size_t n = 0; n < N; n++) {
98  // pull out values of arguments
99  const double y_dbl = value_of(y_vec[n]);
100  const double mu_dbl = value_of(mu_vec[n]);
101 
102  // reusable subexpression values
103  const double y_minus_mu_over_sigma
104  = (y_dbl - mu_dbl) * inv_sigma[n];
105  const double y_minus_mu_over_sigma_squared
106  = y_minus_mu_over_sigma * y_minus_mu_over_sigma;
107 
108  static double NEGATIVE_HALF = - 0.5;
109 
110  // log probability
111  if (include_summand<propto>::value)
112  logp += NEG_LOG_SQRT_TWO_PI;
113  if (include_summand<propto,T_scale>::value)
114  logp -= log_sigma[n];
115  if (include_summand<propto,T_y,T_loc,T_scale>::value)
116  logp += NEGATIVE_HALF * y_minus_mu_over_sigma_squared;
117 
118  // gradients
119  double scaled_diff = inv_sigma[n] * y_minus_mu_over_sigma;
120  if (!is_constant_struct<T_y>::value)
121  operands_and_partials.d_x1[n] -= scaled_diff;
122  if (!is_constant_struct<T_loc>::value)
123  operands_and_partials.d_x2[n] += scaled_diff;
124  if (!is_constant_struct<T_scale>::value)
125  operands_and_partials.d_x3[n]
126  += -inv_sigma[n] + inv_sigma[n] * y_minus_mu_over_sigma_squared;
127  }
128  return operands_and_partials.to_var(logp);
129  }
130 
131 
132  template <bool propto,
133  typename T_y, typename T_loc, typename T_scale>
134  inline
135  typename return_type<T_y,T_loc,T_scale>::type
136  normal_log(const T_y& y, const T_loc& mu, const T_scale& sigma) {
137  return normal_log<propto>(y,mu,sigma,stan::math::default_policy());
138  }
139 
140  template <typename T_y, typename T_loc, typename T_scale,
141  class Policy>
142  inline
143  typename return_type<T_y,T_loc,T_scale>::type
144  normal_log(const T_y& y, const T_loc& mu, const T_scale& sigma,
145  const Policy&) {
146  return normal_log<false>(y,mu,sigma,Policy());
147  }
148 
149  template <typename T_y, typename T_loc, typename T_scale>
150  inline
151  typename return_type<T_y,T_loc,T_scale>::type
152  normal_log(const T_y& y, const T_loc& mu, const T_scale& sigma) {
153  return normal_log<false>(y,mu,sigma,stan::math::default_policy());
154  }
155 
156 
175  template <typename T_y, typename T_loc, typename T_scale,
176  class Policy>
177  typename return_type<T_y,T_loc,T_scale>::type
178  normal_cdf(const T_y& y, const T_loc& mu, const T_scale& sigma,
179  const Policy&) {
180  static const char* function = "stan::prob::normal_cdf(%1%)";
181 
182  using stan::math::check_positive;
183  using stan::math::check_finite;
184  using stan::math::check_not_nan;
185  using stan::math::check_consistent_sizes;
186 
187  typename return_type<T_y, T_loc, T_scale>::type cdf(1);
188  if (!check_not_nan(function, y, "Random variable", &cdf, Policy()))
189  return cdf;
190  if (!check_finite(function, mu, "Location parameter", &cdf, Policy()))
191  return cdf;
192  if (!check_not_nan(function, sigma, "Scale parameter",
193  &cdf, Policy()))
194  return cdf;
195  if (!check_positive(function, sigma, "Scale parameter",
196  &cdf, Policy()))
197  return cdf;
198  if (!(check_consistent_sizes(function,
199  y,mu,sigma,
200  "Random variable","Location parameter","Scale parameter",
201  &cdf, Policy())))
202  return cdf;
203 
204  // check if any vectors are zero length
205  if (!(stan::length(y)
206  && stan::length(mu)
207  && stan::length(sigma)))
208  return 0.0;
209 
210  VectorView<const T_y> y_vec(y);
211  VectorView<const T_loc> mu_vec(mu);
212  VectorView<const T_scale> sigma_vec(sigma);
213  size_t N = max_size(y, mu, sigma);
214 
215  for (size_t n = 0; n < N; n++) {
216  cdf *= 0.5 + 0.5 * erf((y_vec[n] - mu_vec[n]) / (sigma_vec[n] * SQRT_2));
217  }
218  return cdf;
219  }
220 
221  template <typename T_y, typename T_loc, typename T_scale>
222  inline
223  typename return_type<T_y, T_loc, T_scale>::type
224  normal_cdf(const T_y& y, const T_loc& mu, const T_scale& sigma) {
225  return normal_cdf(y,mu,sigma,stan::math::default_policy());
226  }
227 
228 
229  template <typename T_loc, typename T_scale, class RNG>
230  inline double
231  normal_random(const T_loc& mu, const T_scale& sigma, RNG& rng) {
232  using boost::variate_generator;
233  using boost::normal_distribution;
234  using stan::math::value_of;
235  variate_generator<RNG&, normal_distribution<> >
236  rng_unit_norm(rng, normal_distribution<>());
237  return value_of(mu) + value_of(sigma) * rng_unit_norm();
238  }
239  }
240 }
241 #endif
     [ Stan Home Page ] © 2011–2012, Stan Development Team.