Stan  1.1
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 <stan/agrad.hpp>
5 #include <stan/math/error_handling.hpp>
6 #include <stan/math/special_functions.hpp>
7 #include <stan/meta/traits.hpp>
8 #include <stan/prob/constants.hpp>
9 #include <stan/prob/traits.hpp>
10 
11 namespace stan {
12  namespace prob {
13 
14  // Logistic(y|mu,sigma) [sigma > 0]
15  // FIXME: document
16  template <bool propto,
17  typename T_y, typename T_loc, typename T_scale,
18  class Policy>
19  typename return_type<T_y,T_loc,T_scale>::type
20  logistic_log(const T_y& y, const T_loc& mu, const T_scale& sigma,
21  const Policy&) {
22  static const char* function = "stan::prob::logistic_log(%1%)";
23 
24  using stan::math::check_positive;
25  using stan::math::check_finite;
26  using stan::math::check_consistent_sizes;
27  using stan::math::value_of;
28  using stan::prob::include_summand;
29 
30  // check if any vectors are zero length
31  if (!(stan::length(y)
32  && stan::length(mu)
33  && stan::length(sigma)))
34  return 0.0;
35 
36 
37  // set up return value accumulator
38  double logp(0.0);
39 
40  // validate args (here done over var, which should be OK)
41  if (!check_finite(function, y, "Random variable", &logp, Policy()))
42  return logp;
43  if (!check_finite(function, mu, "Location parameter",
44  &logp, Policy()))
45  return logp;
46  if (!check_finite(function, sigma, "Scale parameter", &logp,
47  Policy()))
48  return logp;
49  if (!check_positive(function, sigma, "Scale parameter",
50  &logp, Policy()))
51  return logp;
52  if (!(check_consistent_sizes(function,
53  y,mu,sigma,
54  "Random variable","Location parameter","Scale parameter",
55  &logp, Policy())))
56  return logp;
57 
58  // check if no variables are involved and prop-to
59  if (!include_summand<propto,T_y,T_loc,T_scale>::value)
60  return 0.0;
61 
62 
63  // set up template expressions wrapping scalars into vector views
64  agrad::OperandsAndPartials<T_y, T_loc, T_scale> operands_and_partials(y, mu, sigma);
65 
66  VectorView<const T_y> y_vec(y);
67  VectorView<const T_loc> mu_vec(mu);
68  VectorView<const T_scale> sigma_vec(sigma);
69  size_t N = max_size(y, mu, sigma);
70 
71  DoubleVectorView<true,is_vector<T_scale>::value> inv_sigma(length(sigma));
72  DoubleVectorView<include_summand<propto,T_scale>::value,is_vector<T_scale>::value> log_sigma(length(sigma));
73  for (size_t i = 0; i < length(sigma); i++) {
74  inv_sigma[i] = 1.0 / value_of(sigma_vec[i]);
75  if (include_summand<propto,T_scale>::value)
76  log_sigma[i] = log(value_of(sigma_vec[i]));
77  }
78 
79  DoubleVectorView<!is_constant_struct<T_loc>::value, is_vector<T_loc>::value || is_vector<T_scale>::value>
80  exp_mu_div_sigma(max_size(mu,sigma));
81  DoubleVectorView<!is_constant_struct<T_loc>::value, is_vector<T_y>::value || is_vector<T_scale>::value>
82  exp_y_div_sigma(max_size(y,sigma));
83  if (!is_constant_struct<T_loc>::value) {
84  for (size_t n = 0; n < max_size(mu,sigma); n++)
85  exp_mu_div_sigma[n] = exp(value_of(mu_vec[n]) / value_of(sigma_vec[n]));
86  for (size_t n = 0; n < max_size(y,sigma); n++)
87  exp_y_div_sigma[n] = exp(value_of(y_vec[n]) / value_of(sigma_vec[n]));
88  }
89 
90  using stan::math::log1p;
91  for (size_t n = 0; n < N; n++) {
92  const double y_dbl = value_of(y_vec[n]);
93  const double mu_dbl = value_of(mu_vec[n]);
94 
95  const double y_minus_mu = y_dbl - mu_dbl;
96  const double y_minus_mu_div_sigma = y_minus_mu * inv_sigma[n];
97  double exp_m_y_minus_mu_div_sigma(0);
98  if (include_summand<propto,T_y,T_loc,T_scale>::value)
99  exp_m_y_minus_mu_div_sigma = exp(-y_minus_mu_div_sigma);
100  double inv_1p_exp_y_minus_mu_div_sigma(0);
101  if (!is_constant_struct<T_y>::value || !is_constant_struct<T_scale>::value)
102  inv_1p_exp_y_minus_mu_div_sigma = 1 / (1 + exp(y_minus_mu_div_sigma));
103 
104  if (include_summand<propto,T_y,T_loc,T_scale>::value)
105  logp -= y_minus_mu_div_sigma;
106  if (include_summand<propto,T_scale>::value)
107  logp -= log_sigma[n];
108  if (include_summand<propto,T_y,T_loc,T_scale>::value)
109  logp -= 2.0 * log1p(exp_m_y_minus_mu_div_sigma);
110 
111  if (!is_constant_struct<T_y>::value)
112  operands_and_partials.d_x1[n] += (2 * inv_1p_exp_y_minus_mu_div_sigma - 1) * inv_sigma[n];
113  if (!is_constant_struct<T_loc>::value)
114  operands_and_partials.d_x2[n] +=
115  (1 - 2 * exp_mu_div_sigma[n] / (exp_mu_div_sigma[n] + exp_y_div_sigma[n])) * inv_sigma[n];
116  if (!is_constant_struct<T_scale>::value)
117  operands_and_partials.d_x3[n] +=
118  ((1 - 2 * inv_1p_exp_y_minus_mu_div_sigma)*y_minus_mu*inv_sigma[n] - 1) * inv_sigma[n];
119  }
120  return operands_and_partials.to_var(logp);
121  }
122 
123  template <bool propto,
124  typename T_y, typename T_loc, typename T_scale>
125  inline
126  typename return_type<T_y,T_loc,T_scale>::type
127  logistic_log(const T_y& y, const T_loc& mu, const T_scale& sigma) {
128  return logistic_log<propto>(y,mu,sigma,stan::math::default_policy());
129  }
130 
131  template <typename T_y, typename T_loc, typename T_scale,
132  class Policy>
133  inline
134  typename return_type<T_y,T_loc,T_scale>::type
135  logistic_log(const T_y& y, const T_loc& mu, const T_scale& sigma,
136  const Policy&) {
137  return logistic_log<false>(y,mu,sigma,Policy());
138  }
139 
140  template <typename T_y, typename T_loc, typename T_scale>
141  inline
142  typename return_type<T_y,T_loc,T_scale>::type
143  logistic_log(const T_y& y, const T_loc& mu, const T_scale& sigma) {
144  return logistic_log<false>(y,mu,sigma,stan::math::default_policy());
145  }
146 
147  // Logistic(y|mu,sigma) [sigma > 0]
148  template <typename T_y, typename T_loc, typename T_scale, class Policy>
149  typename return_type<T_y, T_loc, T_scale>::type
150  logistic_cdf(const T_y& y, const T_loc& mu, const T_scale& sigma, const Policy&) {
151 
152  // Size checks
153  if ( !( stan::length(y) && stan::length(mu) && stan::length(sigma) ) ) return 0.0;
154 
155  // Error checks
156  static const char* function = "stan::prob::logistic_cdf(%1%)";
157 
158  using stan::math::check_positive;
159  using stan::math::check_finite;
160  using stan::math::check_consistent_sizes;
161 
162  using boost::math::tools::promote_args;
163 
164  double P(1.0);
165 
166  if (!check_finite(function, y, "Random variable", &P, Policy()))
167  return P;
168 
169  if (!check_finite(function, mu, "Location parameter", &P, Policy()))
170  return P;
171 
172  if (!check_finite(function, sigma, "Scale parameter", &P, Policy()))
173  return P;
174 
175  if (!check_positive(function, sigma, "Scale parameter", &P, Policy()))
176  return P;
177 
178  if (!(check_consistent_sizes(function, y, mu, sigma,
179  "Random variable", "Location parameter", "Scale parameter",
180  &P, Policy())))
181  return P;
182 
183  // Wrap arguments in vectors
184  VectorView<const T_y> y_vec(y);
185  VectorView<const T_loc> mu_vec(mu);
186  VectorView<const T_scale> sigma_vec(sigma);
187  size_t N = max_size(y, mu, sigma);
188 
189  agrad::OperandsAndPartials<T_y, T_loc, T_scale> operands_and_partials(y, mu, sigma);
190 
191  std::fill(operands_and_partials.all_partials,
192  operands_and_partials.all_partials + operands_and_partials.nvaris, 0.0);
193 
194  // Compute vectorized CDF and its gradients
195  using stan::math::value_of;
196 
197  for (size_t n = 0; n < N; n++) {
198 
199  // Pull out values
200  const double y_dbl = value_of(y_vec[n]);
201  const double mu_dbl = value_of(mu_vec[n]);
202  const double sigma_dbl = value_of(sigma_vec[n]);
203  const double sigma_inv_vec = 1.0 / value_of(sigma_vec[n]);
204 
205  // Compute
206  const double Pn = 1.0 / ( 1.0 + exp( - (y_dbl - mu_dbl) * sigma_inv_vec ) );
207 
208  P *= Pn;
209 
210  if (!is_constant_struct<T_y>::value)
211  operands_and_partials.d_x1[n]
212  += exp(logistic_log(y_dbl, mu_dbl, sigma_dbl)) / Pn;
213 
214  if (!is_constant_struct<T_loc>::value)
215  operands_and_partials.d_x2[n]
216  += - exp(logistic_log(y_dbl, mu_dbl, sigma_dbl)) / Pn;
217 
218  if (!is_constant_struct<T_scale>::value)
219  operands_and_partials.d_x3[n]
220  += - (y_dbl - mu_dbl) * sigma_inv_vec * exp(logistic_log(y_dbl, mu_dbl, sigma_dbl)) / Pn;
221 
222  }
223 
224  for (size_t n = 0; n < N; n++) {
225 
226  if (!is_constant_struct<T_y>::value)
227  operands_and_partials.d_x1[n] *= P;
228 
229  if (!is_constant_struct<T_loc>::value)
230  operands_and_partials.d_x2[n] *= P;
231 
232  if (!is_constant_struct<T_scale>::value)
233  operands_and_partials.d_x3[n] *= P;
234 
235  }
236 
237  return operands_and_partials.to_var(P);
238 
239  }
240 
241  template <typename T_y, typename T_loc, typename T_scale>
242  inline typename return_type<T_y, T_loc, T_scale>::type
243  logistic_cdf(const T_y& y, const T_loc& mu, const T_scale& sigma) {
244  return logistic_cdf(y, mu, sigma, stan::math::default_policy());
245  }
246 
247  }
248 }
249 #endif
     [ Stan Home Page ] © 2011–2013, Stan Development Team.