Stan  1.1
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
student_t.hpp
Go to the documentation of this file.
1 #ifndef __STAN__PROB__DISTRIBUTIONS__UNIVARIATE__CONTINUOUS__STUDENT_T_HPP__
2 #define __STAN__PROB__DISTRIBUTIONS__UNIVARIATE__CONTINUOUS__STUDENT_T_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 #include <stan/prob/internal_math.hpp>
11 
12 namespace stan {
13 
14  namespace prob {
15 
41  template <bool propto, typename T_y, typename T_dof,
42  typename T_loc, typename T_scale,
43  class Policy>
44  typename return_type<T_y,T_dof,T_loc,T_scale>::type
45  student_t_log(const T_y& y, const T_dof& nu, const T_loc& mu,
46  const T_scale& sigma,
47  const Policy&) {
48  static const char* function = "stan::prob::student_t_log(%1%)";
49 
50  using stan::math::check_positive;
51  using stan::math::check_finite;
52  using stan::math::check_not_nan;
53  using stan::math::check_consistent_sizes;
54 
55  // check if any vectors are zero length
56  if (!(stan::length(y)
57  && stan::length(nu)
58  && stan::length(mu)
59  && stan::length(sigma)))
60  return 0.0;
61 
62  typename return_type<T_y,T_dof,T_loc,T_scale>::type logp = 0.0;
63 
64  // validate args (here done over var, which should be OK)
65  if (!check_not_nan(function, y, "Random variable", &logp, Policy()))
66  return logp;
67  if(!check_finite(function, nu, "Degrees of freedom parameter", &logp, Policy()))
68  return logp;
69  if(!check_positive(function, nu, "Degrees of freedom parameter", &logp, Policy()))
70  return logp;
71  if (!check_finite(function, mu, "Location parameter",
72  &logp, Policy()))
73  return logp;
74  if (!check_finite(function, sigma, "Scale parameter",
75  &logp, Policy()))
76  return logp;
77  if (!check_positive(function, sigma, "Scale parameter",
78  &logp, Policy()))
79  return logp;
80 
81 
82  if (!(check_consistent_sizes(function,
83  y,nu,mu,sigma,
84  "Random variable","Degrees of freedom parameter","Location parameter","Scale parameter",
85  &logp, Policy())))
86  return logp;
87 
88  // check if no variables are involved and prop-to
89  if (!include_summand<propto,T_y,T_dof,T_loc,T_scale>::value)
90  return 0.0;
91 
92  VectorView<const T_y> y_vec(y);
93  VectorView<const T_dof> nu_vec(nu);
94  VectorView<const T_loc> mu_vec(mu);
95  VectorView<const T_scale> sigma_vec(sigma);
96  size_t N = max_size(y, nu, mu, sigma);
97 
98  using stan::math::square;
99  using boost::math::lgamma;
100 
101  using std::log;
102 
103  for (size_t n = 0; n < N; n++) {
104  if (include_summand<propto,T_dof>::value)
105  logp += lgamma( (nu_vec[n] + 1.0) / 2.0) - lgamma(nu_vec[n] / 2.0);
106  if (include_summand<propto>::value)
107  logp += NEG_LOG_SQRT_PI;
108  if (include_summand<propto,T_dof>::value)
109  logp -= 0.5 * log(nu_vec[n]);
110  if (include_summand<propto,T_scale>::value)
111  logp -= log(sigma_vec[n]);
112  if (include_summand<propto,T_y,T_dof,T_loc,T_scale>::value)
113  logp -= ((nu_vec[n] + 1.0) / 2.0)
114  * log1p( square(((y_vec[n] - mu_vec[n]) / sigma_vec[n])) / nu_vec[n]);
115  }
116  return logp;
117  }
118 
119  template <bool propto,
120  typename T_y, typename T_dof, typename T_loc, typename T_scale>
121  inline
122  typename return_type<T_y,T_dof,T_loc,T_scale>::type
123  student_t_log(const T_y& y, const T_dof& nu, const T_loc& mu,
124  const T_scale& sigma) {
125  return student_t_log<propto>(y,nu,mu,sigma,stan::math::default_policy());
126  }
127 
128  template <typename T_y, typename T_dof, typename T_loc, typename T_scale,
129  class Policy>
130  inline
131  typename return_type<T_y,T_dof,T_loc,T_scale>::type
132  student_t_log(const T_y& y, const T_dof& nu, const T_loc& mu,
133  const T_scale& sigma,
134  const Policy&) {
135  return student_t_log<false>(y,nu,mu,sigma,Policy());
136  }
137 
138  template <typename T_y, typename T_dof, typename T_loc, typename T_scale>
139  inline
140  typename return_type<T_y,T_dof,T_loc,T_scale>::type
141  student_t_log(const T_y& y, const T_dof& nu, const T_loc& mu,
142  const T_scale& sigma) {
143  return student_t_log<false>(y,nu,mu,sigma,stan::math::default_policy());
144  }
145 
146  template <typename T_y, typename T_dof, typename T_loc, typename T_scale, class Policy>
147  typename boost::math::tools::promote_args<T_y, T_dof, T_loc, T_scale>::type
148  student_t_cdf(const T_y& y, const T_dof& nu, const T_loc& mu, const T_scale& sigma, const Policy&) {
149 
150  static const char* function = "stan::prob::student_t_cdf(%1%)";
151 
152  using stan::math::check_positive;
153  using stan::math::check_finite;
154  using stan::math::check_not_nan;
155  using stan::math::check_consistent_sizes;
156 
157  using boost::math::tools::promote_args;
158 
159 
160  double P(1.0);
161 
162  if (!check_not_nan(function, y, "Random variable", &P, Policy()))
163  return P;
164 
165  if(!check_finite(function, nu, "Degrees of freedom parameter", &P, Policy()))
166  return P;
167 
168  if(!check_positive(function, nu, "Degrees of freedom parameter", &P, Policy()))
169  return P;
170 
171  if (!check_finite(function, mu, "Location parameter",
172  &P, Policy()))
173  return P;
174 
175  if (!check_finite(function, sigma, "Scale parameter",
176  &P, Policy()))
177  return P;
178  if (!check_positive(function, sigma, "Scale parameter",
179  &P, Policy()))
180  return P;
181 
182  // Wrap arguments in vectors
183  VectorView<const T_y> y_vec(y);
184  VectorView<const T_dof> nu_vec(nu);
185  VectorView<const T_loc> mu_vec(mu);
186  VectorView<const T_scale> sigma_vec(sigma);
187  //size_t N = max_size(y, nu, mu, sigma);
188  size_t N = max_size(y, nu);
189 
190  agrad::OperandsAndPartials<T_y, T_dof, T_loc, T_scale> operands_and_partials(y, nu, mu, sigma);
191 
192  std::fill(operands_and_partials.all_partials,
193  operands_and_partials.all_partials + operands_and_partials.nvaris, 0.0);
194 
195  using stan::math::value_of;
196  using boost::math::ibeta;
197  using boost::math::ibeta_derivative;
198 
199  using boost::math::digamma;
200  using boost::math::beta;
201 
202  // Cache a few expensive function calls if nu is a parameter
203  double digammaHalf = 0;
204 
205  DoubleVectorView<!is_constant_struct<T_dof>::value,is_vector<T_dof>::value> digamma_vec(stan::length(nu));
206  DoubleVectorView<!is_constant_struct<T_dof>::value,is_vector<T_dof>::value> digammaNu_vec(stan::length(nu));
207  DoubleVectorView<!is_constant_struct<T_dof>::value,is_vector<T_dof>::value> digammaNuPlusHalf_vec(stan::length(nu));
208  DoubleVectorView<!is_constant_struct<T_dof>::value,is_vector<T_dof>::value> betaNuHalf_vec(stan::length(nu));
209 
210  if (!is_constant_struct<T_dof>::value) {
211 
212  digammaHalf = digamma(0.5);
213 
214  for (size_t i = 0; i < stan::length(nu); i++) {
215 
216  const double nu_dbl = value_of(nu_vec[i]);
217 
218  digammaNu_vec[i] = digamma(0.5 * nu_dbl);
219  digammaNuPlusHalf_vec[i] = digamma(0.5 + 0.5 * nu_dbl);
220  betaNuHalf_vec[i] = beta(0.5, 0.5 * nu_dbl);
221 
222  }
223 
224  }
225 
226  // Compute vectorized CDF and gradient
227  for (size_t n = 0; n < N; n++) {
228 
229  const double sigma_inv = 1.0 / value_of(sigma_vec[n]);
230  const double t = (value_of(y_vec[n]) - value_of(mu_vec[n])) * sigma_inv;
231  const double nu_dbl = value_of(nu_vec[n]);
232  const double q = nu_dbl / (t * t);
233  const double r = 1.0 / (1.0 + q);
234  double J = - 2 * r * r * q / t;
235 
236  if(q < 2)
237  {
238 
239  double z = ibeta(0.5 * nu_dbl, 0.5, 1.0 - r);
240 
241  const double Pn = t > 0 ? 1.0 - 0.5 * z : 0.5 * z;
242 
243  double d_ibeta = ibeta_derivative(0.5 * nu_dbl, 0.5, 1.0 - r);
244 
245  P *= Pn;
246 
247  if (!is_constant_struct<T_y>::value)
248  operands_and_partials.d_x1[n]
249  += (d_ibeta * J * sigma_inv) / Pn;
250 
251  if (!is_constant_struct<T_dof>::value) {
252 
253  double g1 = 0;
254  double g2 = 0;
255 
256  stan::math::gradRegIncBeta(g1, g2, 0.5 * nu_dbl, 0.5, 1.0 - r,
257  digammaNu_vec[n], digammaHalf, digammaNuPlusHalf_vec[n], betaNuHalf_vec[n]);
258 
259  operands_and_partials.d_x2[n]
260  += 0.5 * ( d_ibeta * (r / t) * (r / t) + 0.5 * g1 ) / Pn;
261 
262  }
263 
264  if (!is_constant_struct<T_loc>::value)
265  operands_and_partials.d_x3[n]
266  += (- d_ibeta * J * sigma_inv) / Pn;
267 
268  if (!is_constant_struct<T_scale>::value)
269  operands_and_partials.d_x4[n]
270  += (- d_ibeta * J * sigma_inv * t) / Pn;
271 
272  }
273  else
274  {
275 
276  double z = 1.0 - ibeta(0.5, 0.5 * nu_dbl, r);
277 
278  const double Pn = t > 0 ? 1.0 - 0.5 * z : 0.5 * z;
279 
280  double d_ibeta = ibeta_derivative(0.5, 0.5 * nu_dbl, r);
281 
282  P *= Pn;
283 
284  if (!is_constant_struct<T_y>::value)
285  operands_and_partials.d_x1[n]
286  += (d_ibeta * J * sigma_inv) / Pn;
287 
288  if (!is_constant_struct<T_dof>::value) {
289 
290  double g1 = 0;
291  double g2 = 0;
292 
293  stan::math::gradRegIncBeta(g1, g2, 0.5, 0.5 * nu_dbl, r,
294  digammaHalf, digammaNu_vec[n], digammaNuPlusHalf_vec[n], betaNuHalf_vec[n]);
295 
296  operands_and_partials.d_x2[n]
297  += 0.5 * ( - d_ibeta * (r / t) * (r / t) + 0.5 * g2 ) / Pn;
298 
299  }
300 
301  if (!is_constant_struct<T_loc>::value)
302  operands_and_partials.d_x3[n]
303  += (- d_ibeta * J * sigma_inv) / Pn;
304 
305  if (!is_constant_struct<T_scale>::value)
306  operands_and_partials.d_x4[n]
307  += (- d_ibeta * J * sigma_inv * t) / Pn;
308 
309  }
310 
311  if(t > 0) operands_and_partials.d_x2[n] *= -1;
312 
313  }
314 
315  for (size_t n = 0; n < N; n++) {
316 
317  if (!is_constant_struct<T_y>::value)
318  operands_and_partials.d_x1[n] *= P;
319 
320  if (!is_constant_struct<T_dof>::value)
321  operands_and_partials.d_x2[n] *= P;
322 
323  if (!is_constant_struct<T_loc>::value)
324  operands_and_partials.d_x3[n] *= P;
325 
326  if (!is_constant_struct<T_scale>::value)
327  operands_and_partials.d_x4[n] *= P;
328 
329  }
330 
331  return operands_and_partials.to_var(P);
332 
333  }
334 
335  template <typename T_y, typename T_dof, typename T_loc, typename T_scale>
336  typename boost::math::tools::promote_args<T_y,T_dof,T_loc,T_scale>::type
337  student_t_cdf(const T_y& y, const T_dof& nu, const T_loc& mu, const T_scale& sigma) {
338  return student_t_cdf(y, nu, mu, sigma, stan::math::default_policy());
339  }
340 
341  }
342 }
343 #endif
     [ Stan Home Page ] © 2011–2013, Stan Development Team.