Stan  1.3
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
logistic.hpp
Go to the documentation of this file.
1 #ifndef __STAN__PROB__DISTRIBUTIONS__UNIVARIATE__CONTINUOUS__LOGISTIC_HPP__
2 #define __STAN__PROB__DISTRIBUTIONS__UNIVARIATE__CONTINUOUS__LOGISTIC_HPP__
3 
4 #include <boost/random/exponential_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/meta/traits.hpp>
11 #include <stan/prob/constants.hpp>
12 #include <stan/prob/traits.hpp>
13 
14 namespace stan {
15  namespace prob {
16 
17  // Logistic(y|mu,sigma) [sigma > 0]
18  // FIXME: document
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  logistic_log(const T_y& y, const T_loc& mu, const T_scale& sigma,
24  const Policy&) {
25  static const char* function = "stan::prob::logistic_log(%1%)";
26 
27  using stan::math::check_positive;
28  using stan::math::check_finite;
29  using stan::math::check_consistent_sizes;
30  using stan::math::value_of;
31  using stan::prob::include_summand;
32 
33  // check if any vectors are zero length
34  if (!(stan::length(y)
35  && stan::length(mu)
36  && stan::length(sigma)))
37  return 0.0;
38 
39 
40  // set up return value accumulator
41  double logp(0.0);
42 
43  // validate args (here done over var, which should be OK)
44  if (!check_finite(function, y, "Random variable", &logp, Policy()))
45  return logp;
46  if (!check_finite(function, mu, "Location parameter",
47  &logp, Policy()))
48  return logp;
49  if (!check_finite(function, sigma, "Scale parameter", &logp,
50  Policy()))
51  return logp;
52  if (!check_positive(function, sigma, "Scale parameter",
53  &logp, Policy()))
54  return logp;
55  if (!(check_consistent_sizes(function,
56  y,mu,sigma,
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 
66  // set up template expressions wrapping scalars into vector views
67  agrad::OperandsAndPartials<T_y, T_loc, T_scale> operands_and_partials(y, mu, sigma);
68 
69  VectorView<const T_y> y_vec(y);
70  VectorView<const T_loc> mu_vec(mu);
71  VectorView<const T_scale> sigma_vec(sigma);
72  size_t N = max_size(y, mu, sigma);
73 
74  DoubleVectorView<true,is_vector<T_scale>::value> inv_sigma(length(sigma));
75  DoubleVectorView<include_summand<propto,T_scale>::value,is_vector<T_scale>::value> log_sigma(length(sigma));
76  for (size_t i = 0; i < length(sigma); i++) {
77  inv_sigma[i] = 1.0 / value_of(sigma_vec[i]);
78  if (include_summand<propto,T_scale>::value)
79  log_sigma[i] = log(value_of(sigma_vec[i]));
80  }
81 
82  DoubleVectorView<!is_constant_struct<T_loc>::value, is_vector<T_loc>::value || is_vector<T_scale>::value>
83  exp_mu_div_sigma(max_size(mu,sigma));
84  DoubleVectorView<!is_constant_struct<T_loc>::value, is_vector<T_y>::value || is_vector<T_scale>::value>
85  exp_y_div_sigma(max_size(y,sigma));
86  if (!is_constant_struct<T_loc>::value) {
87  for (size_t n = 0; n < max_size(mu,sigma); n++)
88  exp_mu_div_sigma[n] = exp(value_of(mu_vec[n]) / value_of(sigma_vec[n]));
89  for (size_t n = 0; n < max_size(y,sigma); n++)
90  exp_y_div_sigma[n] = exp(value_of(y_vec[n]) / value_of(sigma_vec[n]));
91  }
92 
93  using stan::math::log1p;
94  for (size_t n = 0; n < N; n++) {
95  const double y_dbl = value_of(y_vec[n]);
96  const double mu_dbl = value_of(mu_vec[n]);
97 
98  const double y_minus_mu = y_dbl - mu_dbl;
99  const double y_minus_mu_div_sigma = y_minus_mu * inv_sigma[n];
100  double exp_m_y_minus_mu_div_sigma(0);
101  if (include_summand<propto,T_y,T_loc,T_scale>::value)
102  exp_m_y_minus_mu_div_sigma = exp(-y_minus_mu_div_sigma);
103  double inv_1p_exp_y_minus_mu_div_sigma(0);
104  if (!is_constant_struct<T_y>::value || !is_constant_struct<T_scale>::value)
105  inv_1p_exp_y_minus_mu_div_sigma = 1 / (1 + exp(y_minus_mu_div_sigma));
106 
107  if (include_summand<propto,T_y,T_loc,T_scale>::value)
108  logp -= y_minus_mu_div_sigma;
109  if (include_summand<propto,T_scale>::value)
110  logp -= log_sigma[n];
111  if (include_summand<propto,T_y,T_loc,T_scale>::value)
112  logp -= 2.0 * log1p(exp_m_y_minus_mu_div_sigma);
113 
114  if (!is_constant_struct<T_y>::value)
115  operands_and_partials.d_x1[n] += (2 * inv_1p_exp_y_minus_mu_div_sigma - 1) * inv_sigma[n];
116  if (!is_constant_struct<T_loc>::value)
117  operands_and_partials.d_x2[n] +=
118  (1 - 2 * exp_mu_div_sigma[n] / (exp_mu_div_sigma[n] + exp_y_div_sigma[n])) * inv_sigma[n];
119  if (!is_constant_struct<T_scale>::value)
120  operands_and_partials.d_x3[n] +=
121  ((1 - 2 * inv_1p_exp_y_minus_mu_div_sigma)*y_minus_mu*inv_sigma[n] - 1) * inv_sigma[n];
122  }
123  return operands_and_partials.to_var(logp);
124  }
125 
126  template <bool propto,
127  typename T_y, typename T_loc, typename T_scale>
128  inline
129  typename return_type<T_y,T_loc,T_scale>::type
130  logistic_log(const T_y& y, const T_loc& mu, const T_scale& sigma) {
131  return logistic_log<propto>(y,mu,sigma,stan::math::default_policy());
132  }
133 
134  template <typename T_y, typename T_loc, typename T_scale,
135  class Policy>
136  inline
137  typename return_type<T_y,T_loc,T_scale>::type
138  logistic_log(const T_y& y, const T_loc& mu, const T_scale& sigma,
139  const Policy&) {
140  return logistic_log<false>(y,mu,sigma,Policy());
141  }
142 
143  template <typename T_y, typename T_loc, typename T_scale>
144  inline
145  typename return_type<T_y,T_loc,T_scale>::type
146  logistic_log(const T_y& y, const T_loc& mu, const T_scale& sigma) {
147  return logistic_log<false>(y,mu,sigma,stan::math::default_policy());
148  }
149 
150  // Logistic(y|mu,sigma) [sigma > 0]
151  template <typename T_y, typename T_loc, typename T_scale, class Policy>
152  typename return_type<T_y, T_loc, T_scale>::type
153  logistic_cdf(const T_y& y, const T_loc& mu, const T_scale& sigma, const Policy&) {
154 
155  // Size checks
156  if ( !( stan::length(y) && stan::length(mu) && stan::length(sigma) ) ) return 1.0;
157 
158  // Error checks
159  static const char* function = "stan::prob::logistic_cdf(%1%)";
160 
161  using stan::math::check_not_nan;
162  using stan::math::check_positive;
163  using stan::math::check_finite;
164  using stan::math::check_consistent_sizes;
165  using stan::math::value_of;
166 
167  using boost::math::tools::promote_args;
168 
169  double P(1.0);
170 
171  if (!check_not_nan(function, y, "Random variable", &P, Policy()))
172  return P;
173 
174  if (!check_finite(function, mu, "Location parameter", &P, Policy()))
175  return P;
176 
177  if (!check_finite(function, sigma, "Scale parameter", &P, Policy()))
178  return P;
179 
180  if (!check_positive(function, sigma, "Scale parameter", &P, Policy()))
181  return P;
182 
183  if (!(check_consistent_sizes(function, y, mu, sigma,
184  "Random variable", "Location parameter", "Scale parameter",
185  &P, Policy())))
186  return P;
187 
188  // Wrap arguments in vectors
189  VectorView<const T_y> y_vec(y);
190  VectorView<const T_loc> mu_vec(mu);
191  VectorView<const T_scale> sigma_vec(sigma);
192  size_t N = max_size(y, mu, sigma);
193 
194  agrad::OperandsAndPartials<T_y, T_loc, T_scale> operands_and_partials(y, mu, sigma);
195 
196  std::fill(operands_and_partials.all_partials,
197  operands_and_partials.all_partials + operands_and_partials.nvaris, 0.0);
198 
199  // Explicit return for extreme values
200  // The gradients are technically ill-defined, but treated as zero
201 
202  for (size_t i = 0; i < stan::length(y); i++) {
203  if (value_of(y_vec[i]) == -std::numeric_limits<double>::infinity())
204  return operands_and_partials.to_var(0.0);
205  }
206 
207  // Compute vectorized CDF and its gradients
208  for (size_t n = 0; n < N; n++) {
209 
210  // Explicit results for extreme values
211  // The gradients are technically ill-defined, but treated as zero
212  if (value_of(y_vec[n]) == std::numeric_limits<double>::infinity()) {
213  continue;
214  }
215 
216  // Pull out values
217  const double y_dbl = value_of(y_vec[n]);
218  const double mu_dbl = value_of(mu_vec[n]);
219  const double sigma_dbl = value_of(sigma_vec[n]);
220  const double sigma_inv_vec = 1.0 / value_of(sigma_vec[n]);
221 
222  // Compute
223  const double Pn = 1.0 / ( 1.0 + exp( - (y_dbl - mu_dbl) * sigma_inv_vec ) );
224 
225  P *= Pn;
226 
227  if (!is_constant_struct<T_y>::value)
228  operands_and_partials.d_x1[n]
229  += exp(logistic_log(y_dbl, mu_dbl, sigma_dbl)) / Pn;
230 
231  if (!is_constant_struct<T_loc>::value)
232  operands_and_partials.d_x2[n]
233  += - exp(logistic_log(y_dbl, mu_dbl, sigma_dbl)) / Pn;
234 
235  if (!is_constant_struct<T_scale>::value)
236  operands_and_partials.d_x3[n]
237  += - (y_dbl - mu_dbl) * sigma_inv_vec * exp(logistic_log(y_dbl, mu_dbl, sigma_dbl)) / Pn;
238 
239  }
240 
241  if (!is_constant_struct<T_y>::value) {
242  for(size_t n = 0; n < stan::length(y); ++n) operands_and_partials.d_x1[n] *= P;
243  }
244 
245  if (!is_constant_struct<T_loc>::value) {
246  for(size_t n = 0; n < stan::length(mu); ++n) operands_and_partials.d_x2[n] *= P;
247  }
248 
249  if (!is_constant_struct<T_scale>::value) {
250  for(size_t n = 0; n < stan::length(sigma); ++n) operands_and_partials.d_x3[n] *= P;
251  }
252 
253  return operands_and_partials.to_var(P);
254 
255  }
256 
257  template <typename T_y, typename T_loc, typename T_scale>
258  inline typename return_type<T_y, T_loc, T_scale>::type
259  logistic_cdf(const T_y& y, const T_loc& mu, const T_scale& sigma) {
260  return logistic_cdf(y, mu, sigma, stan::math::default_policy());
261  }
262 
263  template <class RNG>
264  inline double
265  logistic_rng(const double mu,
266  const double sigma,
267  RNG& rng) {
268  using boost::variate_generator;
269  using boost::random::exponential_distribution;
270  variate_generator<RNG&, exponential_distribution<> >
271  exp_rng(rng, exponential_distribution<>(1));
272  return mu - sigma * std::log(exp_rng() / exp_rng());
273  }
274  }
275 }
276 #endif
     [ Stan Home Page ] © 2011–2013, Stan Development Team.