Stan  1.3
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
lognormal.hpp
Go to the documentation of this file.
1 #ifndef __STAN__PROB__DISTRIBUTIONS__UNIVARIATE__CONTINUOUS__LOGNORMAL_HPP__
2 #define __STAN__PROB__DISTRIBUTIONS__UNIVARIATE__CONTINUOUS__LOGNORMAL_HPP__
3 
4 #include <boost/random/lognormal_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/functions/value_of.hpp>
10 #include <stan/math/functions/square.hpp>
11 #include <stan/meta/traits.hpp>
12 #include <stan/prob/constants.hpp>
13 #include <stan/prob/traits.hpp>
14 
15 namespace stan {
16  namespace prob {
17 
18  // LogNormal(y|mu,sigma) [y >= 0; sigma > 0]
19  // FIXME: document
20  template <bool propto,
21  typename T_y, typename T_loc, typename T_scale,
22  class Policy>
23  typename return_type<T_y,T_loc,T_scale>::type
24  lognormal_log(const T_y& y, const T_loc& mu, const T_scale& sigma,
25  const Policy&) {
26  static const char* function = "stan::prob::lognormal_log(%1%)";
27 
28  using stan::is_constant_struct;
29  using stan::math::check_not_nan;
30  using stan::math::check_finite;
31  using stan::math::check_positive;
32  using stan::math::check_consistent_sizes;
33  using stan::math::value_of;
34  using stan::prob::include_summand;
35 
36 
37  // check if any vectors are zero length
38  if (!(stan::length(y)
39  && stan::length(mu)
40  && stan::length(sigma)))
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, sigma, "Scale parameter",
53  &logp, Policy()))
54  return logp;
55  if (!check_positive(function, sigma, "Scale parameter",
56  &logp, Policy()))
57  return logp;
58  if (!(check_consistent_sizes(function,
59  y,mu,sigma,
60  "Random variable","Location parameter","Scale parameter",
61  &logp, Policy())))
62  return logp;
63 
64 
65  VectorView<const T_y> y_vec(y);
66  VectorView<const T_loc> mu_vec(mu);
67  VectorView<const T_scale> sigma_vec(sigma);
68  size_t N = max_size(y, mu, sigma);
69 
70  for (size_t n = 0; n < length(y); n++)
71  if (value_of(y_vec[n]) <= 0)
72  return LOG_ZERO;
73 
74  agrad::OperandsAndPartials<T_y, T_loc, T_scale> operands_and_partials(y, mu, sigma);
75 
76  using stan::math::square;
77  using std::log;
78  using stan::prob::NEG_LOG_SQRT_TWO_PI;
79 
80 
81  DoubleVectorView<include_summand<propto,T_scale>::value,is_vector<T_scale>::value> log_sigma(length(sigma));
82  if (include_summand<propto, T_scale>::value)
83  for (size_t n = 0; n < length(sigma); n++)
84  log_sigma[n] = log(value_of(sigma_vec[n]));
85  DoubleVectorView<include_summand<propto,T_y,T_loc,T_scale>::value,is_vector<T_scale>::value> inv_sigma(length(sigma));
86  DoubleVectorView<include_summand<propto,T_y,T_loc,T_scale>::value,is_vector<T_scale>::value> inv_sigma_sq(length(sigma));
87  if (include_summand<propto,T_y,T_loc,T_scale>::value)
88  for (size_t n = 0; n < length(sigma); n++)
89  inv_sigma[n] = 1 / value_of(sigma_vec[n]);
90  if (include_summand<propto,T_y,T_loc,T_scale>::value)
91  for (size_t n = 0; n < length(sigma); n++)
92  inv_sigma_sq[n] = inv_sigma[n] * inv_sigma[n];
93 
94  DoubleVectorView<include_summand<propto,T_y,T_loc,T_scale>::value,is_vector<T_y>::value> log_y(length(y));
95  if (include_summand<propto,T_y,T_loc,T_scale>::value)
96  for (size_t n = 0; n < length(y); n++)
97  log_y[n] = log(value_of(y_vec[n]));
98  DoubleVectorView<!is_constant_struct<T_y>::value,is_vector<T_y>::value> inv_y(length(y));
99  if (!is_constant_struct<T_y>::value)
100  for (size_t n = 0; n < length(y); n++)
101  inv_y[n] = 1 / value_of(y_vec[n]);
102 
103  if (include_summand<propto>::value)
104  logp += N * NEG_LOG_SQRT_TWO_PI;
105 
106  for (size_t n = 0; n < N; n++) {
107  const double mu_dbl = value_of(mu_vec[n]);
108 
109  double logy_m_mu(0);
110  if (include_summand<propto,T_y,T_loc,T_scale>::value ||
111  !is_constant_struct<T_y>::value)
112  logy_m_mu = log_y[n] - mu_dbl;
113 
114  double logy_m_mu_sq = logy_m_mu * logy_m_mu;
115  double logy_m_mu_div_sigma(0);
116  if (!is_constant_struct<T_y>::value ||
117  !is_constant_struct<T_loc>::value ||
118  !is_constant_struct<T_scale>::value)
119  logy_m_mu_div_sigma = logy_m_mu * inv_sigma_sq[n];
120 
121 
122  // log probability
123  if (include_summand<propto,T_scale>::value)
124  logp -= log_sigma[n];
125  if (include_summand<propto,T_y>::value)
126  logp -= log_y[n];
127  if (include_summand<propto,T_y,T_loc,T_scale>::value)
128  logp -= 0.5 * logy_m_mu_sq * inv_sigma_sq[n];
129 
130  // gradients
131  if (!is_constant_struct<T_y>::value)
132  operands_and_partials.d_x1[n] -= (1 + logy_m_mu_div_sigma) * inv_y[n];
133  if (!is_constant_struct<T_loc>::value)
134  operands_and_partials.d_x2[n] += logy_m_mu_div_sigma;
135  if (!is_constant_struct<T_scale>::value)
136  operands_and_partials.d_x3[n] += (logy_m_mu_div_sigma * logy_m_mu - 1) * inv_sigma[n];
137  }
138  return operands_and_partials.to_var(logp);
139  }
140 
141  template <bool propto,
142  typename T_y, typename T_loc, typename T_scale>
143  inline
144  typename return_type<T_y,T_loc,T_scale>::type
145  lognormal_log(const T_y& y, const T_loc& mu, const T_scale& sigma) {
146  return lognormal_log<propto>(y,mu,sigma,stan::math::default_policy());
147  }
148 
149  template <typename T_y, typename T_loc, typename T_scale,
150  class Policy>
151  inline
152  typename return_type<T_y,T_loc,T_scale>::type
153  lognormal_log(const T_y& y, const T_loc& mu, const T_scale& sigma,
154  const Policy&) {
155  return lognormal_log<false>(y,mu,sigma,Policy());
156  }
157 
158  template <typename T_y, typename T_loc, typename T_scale>
159  inline
160  typename return_type<T_y,T_loc,T_scale>::type
161  lognormal_log(const T_y& y, const T_loc& mu, const T_scale& sigma) {
162  return lognormal_log<false>(y,mu,sigma,stan::math::default_policy());
163  }
164 
165 
166 
167 
168  template <typename T_y, typename T_loc, typename T_scale,
169  class Policy>
170  typename boost::math::tools::promote_args<T_y,T_loc,T_scale>::type
171  lognormal_cdf(const T_y& y, const T_loc& mu, const T_scale& sigma,
172  const Policy&) {
173  static const char* function = "stan::prob::lognormal_cdf(%1%)";
174 
175  using stan::math::check_not_nan;
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;
181  if (!check_not_nan(function, y, "Random variable", &lp, Policy()))
182  return lp;
183  if (!check_finite(function, mu, "Location parameter", &lp, Policy()))
184  return lp;
185  if (!check_finite(function, sigma, "Scale parameter",
186  &lp, Policy()))
187  return lp;
188  if (!check_positive(function, sigma, "Scale parameter",
189  &lp, Policy()))
190  return lp;
191 
192  return 0.5 * erfc(-(log(y) - mu)/(sigma * SQRT_2));
193  }
194 
195  template <typename T_y, typename T_loc, typename T_scale>
196  inline
197  typename boost::math::tools::promote_args<T_y,T_loc,T_scale>::type
198  lognormal_cdf(const T_y& y, const T_loc& mu, const T_scale& sigma) {
199  return lognormal_cdf(y,mu,sigma,stan::math::default_policy());
200  }
201 
202  template <class RNG>
203  inline double
204  lognormal_rng(const double mu,
205  const double sigma,
206  RNG& rng) {
207  using boost::variate_generator;
208  using boost::random::lognormal_distribution;
209  variate_generator<RNG&, lognormal_distribution<> >
210  lognorm_rng(rng, lognormal_distribution<>(mu, sigma));
211  return lognorm_rng();
212  }
213  }
214 }
215 #endif
     [ Stan Home Page ] © 2011–2013, Stan Development Team.