Stan  1.3
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
weibull.hpp
Go to the documentation of this file.
1 #ifndef __STAN__PROB__DISTRIBUTIONS__WEIBULL_HPP__
2 #define __STAN__PROB__DISTRIBUTIONS__WEIBULL_HPP__
3 
4 #include <boost/random/weibull_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 
16  namespace prob {
17 
18  // Weibull(y|sigma,alpha) [y >= 0; sigma > 0; alpha > 0]
19  // FIXME: document
20  template <bool propto,
21  typename T_y, typename T_shape, typename T_scale,
22  class Policy>
23  typename return_type<T_y,T_shape,T_scale>::type
24  weibull_log(const T_y& y, const T_shape& alpha, const T_scale& sigma,
25  const Policy&) {
26  static const char* function = "stan::prob::weibull_log(%1%)";
27 
28  using stan::math::check_finite;
29  using stan::math::check_not_nan;
30  using stan::math::check_positive;
31  using stan::math::value_of;
32  using stan::math::check_consistent_sizes;
33  using stan::math::multiply_log;
34 
35  // check if any vectors are zero length
36  if (!(stan::length(y)
37  && stan::length(alpha)
38  && stan::length(sigma)))
39  return 0.0;
40 
41  // set up return value accumulator
42  double logp(0.0);
43  if(!check_finite(function, y, "Random variable", &logp, Policy()))
44  return logp;
45  if(!check_finite(function, alpha, "Shape parameter",
46  &logp, Policy()))
47  return logp;
48  if(!check_positive(function, alpha, "Shape parameter",
49  &logp, Policy()))
50  return logp;
51  if(!check_not_nan(function, sigma, "Scale parameter",
52  &logp, Policy()))
53  return logp;
54  if(!check_positive(function, sigma, "Scale parameter",
55  &logp, Policy()))
56  return logp;
57  if (!(check_consistent_sizes(function,
58  y,alpha,sigma,
59  "Random variable","Shape parameter","Scale parameter",
60  &logp, Policy())))
61  return logp;
62 
63  // check if no variables are involved and prop-to
64  if (!include_summand<propto,T_y,T_shape,T_scale>::value)
65  return 0.0;
66 
67  VectorView<const T_y> y_vec(y);
68  VectorView<const T_shape> alpha_vec(alpha);
69  VectorView<const T_scale> sigma_vec(sigma);
70  size_t N = max_size(y, alpha, sigma);
71 
72  for (size_t n = 0; n < N; n++) {
73  const double y_dbl = value_of(y_vec[n]);
74  if (y_dbl < 0)
75  return LOG_ZERO;
76  }
77 
78  DoubleVectorView<include_summand<propto,T_shape>::value,
79  is_vector<T_shape>::value> log_alpha(length(alpha));
80  for (size_t i = 0; i < length(alpha); i++)
81  if (include_summand<propto,T_shape>::value)
82  log_alpha[i] = log(value_of(alpha_vec[i]));
83 
84  DoubleVectorView<include_summand<propto,T_y,T_shape>::value,
85  is_vector<T_y>::value> log_y(length(y));
86  for (size_t i = 0; i < length(y); i++)
87  if (include_summand<propto,T_y,T_shape>::value)
88  log_y[i] = log(value_of(y_vec[i]));
89 
90  DoubleVectorView<include_summand<propto,T_shape,T_scale>::value,
91  is_vector<T_scale>::value> log_sigma(length(sigma));
92  for (size_t i = 0; i < length(sigma); i++)
93  if (include_summand<propto,T_shape,T_scale>::value)
94  log_sigma[i] = log(value_of(sigma_vec[i]));
95 
96  DoubleVectorView<include_summand<propto,T_y,T_shape,T_scale>::value,
97  is_vector<T_scale>::value> inv_sigma(length(sigma));
98  for (size_t i = 0; i < length(sigma); i++)
99  if (include_summand<propto,T_y,T_shape,T_scale>::value)
100  inv_sigma[i] = 1.0 / value_of(sigma_vec[i]);
101 
102  DoubleVectorView<include_summand<propto,T_y,T_shape,T_scale>::value,
103  is_vector<T_y>::value | is_vector<T_shape>::value | is_vector<T_scale>::value>
104  y_div_sigma_pow_alpha(N);
105  for (size_t i = 0; i < N; i++)
106  if (include_summand<propto,T_y,T_shape,T_scale>::value) {
107  const double y_dbl = value_of(y_vec[i]);
108  const double alpha_dbl = value_of(alpha_vec[i]);
109  y_div_sigma_pow_alpha[i] = pow(y_dbl * inv_sigma[i], alpha_dbl);
110  }
111 
112  agrad::OperandsAndPartials<T_y,T_shape,T_scale> operands_and_partials(y,alpha,sigma);
113  for (size_t n = 0; n < N; n++) {
114  const double alpha_dbl = value_of(alpha_vec[n]);
115  if (include_summand<propto,T_shape>::value)
116  logp += log_alpha[n];
117  if (include_summand<propto,T_y,T_shape>::value)
118  logp += (alpha_dbl-1.0)*log_y[n];
119  if (include_summand<propto,T_shape,T_scale>::value)
120  logp -= alpha_dbl*log_sigma[n];
121  if (include_summand<propto,T_y,T_shape,T_scale>::value)
122  logp -= y_div_sigma_pow_alpha[n];
123 
124  if (!is_constant_struct<T_y>::value) {
125  const double inv_y = 1.0 / value_of(y_vec[n]);
126  operands_and_partials.d_x1[n]
127  += (alpha_dbl-1.0) * inv_y
128  - alpha_dbl * y_div_sigma_pow_alpha[n] * inv_y;
129  }
130  if (!is_constant_struct<T_shape>::value)
131  operands_and_partials.d_x2[n]
132  += 1.0/alpha_dbl
133  + (1.0 - y_div_sigma_pow_alpha[n]) * (log_y[n] - log_sigma[n]);
134  if (!is_constant_struct<T_scale>::value)
135  operands_and_partials.d_x3[n]
136  += -alpha_dbl * inv_sigma[n]
137  + alpha_dbl * inv_sigma[n] * y_div_sigma_pow_alpha[n];
138  }
139  return operands_and_partials.to_var(logp);
140  }
141 
142 
143  template <bool propto,
144  typename T_y, typename T_shape, typename T_scale>
145  inline
146  typename return_type<T_y,T_shape,T_scale>::type
147  weibull_log(const T_y& y, const T_shape& alpha, const T_scale& sigma) {
148  return weibull_log<propto>(y,alpha,sigma,stan::math::default_policy());
149  }
150 
151 
152  template <typename T_y, typename T_shape, typename T_scale,
153  class Policy>
154  inline
155  typename return_type<T_y,T_shape,T_scale>::type
156  weibull_log(const T_y& y, const T_shape& alpha, const T_scale& sigma,
157  const Policy&) {
158  return weibull_log<false>(y,alpha,sigma,Policy());
159  }
160 
161 
162  template <typename T_y, typename T_shape, typename T_scale>
163  inline
164  typename return_type<T_y,T_shape,T_scale>::type
165  weibull_log(const T_y& y, const T_shape& alpha, const T_scale& sigma) {
166  return weibull_log<false>(y,alpha,sigma,stan::math::default_policy());
167  }
168 
169 
170 
171 
172  template <typename T_y, typename T_shape, typename T_scale,
173  class Policy>
174  typename boost::math::tools::promote_args<T_y,T_shape,T_scale>::type
175  weibull_cdf(const T_y& y, const T_shape& alpha, const T_scale& sigma,
176  const Policy&) {
177 
178  static const char* function = "stan::prob::weibull_cdf(%1%)";
179 
180  using stan::math::check_finite;
181  using stan::math::check_positive;
182  using boost::math::tools::promote_args;
183 
184  typename promote_args<T_y,T_shape,T_scale>::type lp;
185  if (!check_finite(function, alpha, "Shape parameter",
186  &lp, Policy()))
187  return lp;
188  if (!check_positive(function, alpha, "Shape parameter",
189  &lp, Policy()))
190  return lp;
191  if (!check_finite(function, sigma, "Scale parameter",
192  &lp, Policy()))
193  return lp;
194  if (!check_positive(function, sigma, "Scale parameter",
195  &lp, Policy()))
196  return lp;
197 
198  if (y < 0.0)
199  return 0.0;
200  return 1.0 - exp(-pow(y / sigma, alpha));
201  }
202 
203  template <typename T_y, typename T_shape, typename T_scale>
204  inline
205  typename boost::math::tools::promote_args<T_y,T_shape,T_scale>::type
206  weibull_cdf(const T_y& y, const T_shape& alpha, const T_scale& sigma) {
207  return weibull_cdf(y,alpha,sigma,stan::math::default_policy());
208  }
209 
210 
211  template <class RNG>
212  inline double
213  weibull_rng(const double alpha,
214  const double sigma,
215  RNG& rng) {
216  using boost::variate_generator;
217  using boost::random::weibull_distribution;
218  variate_generator<RNG&, weibull_distribution<> >
219  weibull_rng(rng, weibull_distribution<>(alpha, sigma));
220  return weibull_rng();
221  }
222  }
223 }
224 #endif
     [ Stan Home Page ] © 2011–2013, Stan Development Team.