Stan  1.3
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 <boost/random/student_t_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/square.hpp>
10 #include <stan/math/functions/value_of.hpp>
11 #include <stan/meta/traits.hpp>
12 #include <stan/prob/constants.hpp>
13 #include <stan/prob/traits.hpp>
14 #include <stan/prob/internal_math.hpp>
15 
16 namespace stan {
17 
18  namespace prob {
19 
45  template <bool propto, typename T_y, typename T_dof,
46  typename T_loc, typename T_scale,
47  class Policy>
48  typename return_type<T_y,T_dof,T_loc,T_scale>::type
49  student_t_log(const T_y& y, const T_dof& nu, const T_loc& mu,
50  const T_scale& sigma,
51  const Policy&) {
52  static const char* function = "stan::prob::student_t_log(%1%)";
53 
54  using stan::math::check_positive;
55  using stan::math::check_finite;
56  using stan::math::check_not_nan;
57  using stan::math::check_consistent_sizes;
58 
59  // check if any vectors are zero length
60  if (!(stan::length(y)
61  && stan::length(nu)
62  && stan::length(mu)
63  && stan::length(sigma)))
64  return 0.0;
65 
66  double logp(0.0);
67 
68  // validate args (here done over var, which should be OK)
69  if (!check_not_nan(function, y, "Random variable", &logp, Policy()))
70  return logp;
71  if(!check_finite(function, nu, "Degrees of freedom parameter",
72  &logp, Policy()))
73  return logp;
74  if(!check_positive(function, nu, "Degrees of freedom parameter",
75  &logp, Policy()))
76  return logp;
77  if (!check_finite(function, mu, "Location parameter",
78  &logp, Policy()))
79  return logp;
80  if (!check_finite(function, sigma, "Scale parameter",
81  &logp, Policy()))
82  return logp;
83  if (!check_positive(function, sigma, "Scale parameter",
84  &logp, Policy()))
85  return logp;
86 
87 
88  if (!(check_consistent_sizes(function,
89  y,nu,mu,sigma,
90  "Random variable",
91  "Degrees of freedom parameter",
92  "Location parameter","Scale parameter",
93  &logp, Policy())))
94  return logp;
95 
96  // check if no variables are involved and prop-to
97  if (!include_summand<propto,T_y,T_dof,T_loc,T_scale>::value)
98  return 0.0;
99 
100  VectorView<const T_y> y_vec(y);
101  VectorView<const T_dof> nu_vec(nu);
102  VectorView<const T_loc> mu_vec(mu);
103  VectorView<const T_scale> sigma_vec(sigma);
104  size_t N = max_size(y, nu, mu, sigma);
105 
106  using std::log;
107  using boost::math::digamma;
108  using boost::math::lgamma;
109  using stan::math::square;
110  using stan::math::value_of;
111 
112  DoubleVectorView<include_summand<propto,T_y,T_dof,T_loc,T_scale>::value,
113  is_vector<T_dof>::value> half_nu(length(nu));
114  for (size_t i = 0; i < length(nu); i++)
115  if (include_summand<propto,T_y,T_dof,T_loc,T_scale>::value)
116  half_nu[i] = 0.5 * value_of(nu_vec[i]);
117  DoubleVectorView<include_summand<propto,T_dof>::value,
118  is_vector<T_dof>::value> lgamma_half_nu(length(nu));
119  DoubleVectorView<include_summand<propto,T_dof>::value,
120  is_vector<T_dof>::value> lgamma_half_nu_plus_half(length(nu));
121  if (include_summand<propto,T_dof>::value)
122  for (size_t i = 0; i < length(nu); i++) {
123  lgamma_half_nu[i] = lgamma(half_nu[i]);
124  lgamma_half_nu_plus_half[i] = lgamma(half_nu[i] + 0.5);
125  }
126  DoubleVectorView<!is_constant_struct<T_dof>::value,
127  is_vector<T_dof>::value> digamma_half_nu(length(nu));
128  DoubleVectorView<!is_constant_struct<T_dof>::value,
129  is_vector<T_dof>::value> digamma_half_nu_plus_half(length(nu));
130  if (!is_constant_struct<T_dof>::value)
131  for (size_t i = 0; i < length(nu); i++) {
132  digamma_half_nu[i] = digamma(half_nu[i]);
133  digamma_half_nu_plus_half[i] = digamma(half_nu[i] + 0.5);
134  }
135 
136 
137 
138  DoubleVectorView<include_summand<propto,T_dof>::value,
139  is_vector<T_dof>::value> log_nu(length(nu));
140  for (size_t i = 0; i < length(nu); i++)
141  if (include_summand<propto,T_dof>::value)
142  log_nu[i] = log(value_of(nu_vec[i]));
143  DoubleVectorView<include_summand<propto,T_scale>::value,
144  is_vector<T_scale>::value> log_sigma(length(sigma));
145  for (size_t i = 0; i < length(sigma); i++)
146  if (include_summand<propto,T_scale>::value)
147  log_sigma[i] = log(value_of(sigma_vec[i]));
148 
149  DoubleVectorView<include_summand<propto,T_y,T_dof,T_loc,T_scale>::value,
150  is_vector<T_y>::value
151  || is_vector<T_dof>::value
152  || is_vector<T_loc>::value
153  || is_vector<T_scale>::value>
154  square_y_minus_mu_over_sigma__over_nu(N);
155 
156  DoubleVectorView<include_summand<propto,T_y,T_dof,T_loc,T_scale>::value,
157  is_vector<T_y>::value
158  || is_vector<T_dof>::value
159  || is_vector<T_loc>::value
160  || is_vector<T_scale>::value>
161  log1p_exp(N);
162 
163  for (size_t i = 0; i < N; i++)
164  if (include_summand<propto,T_y,T_dof,T_loc,T_scale>::value) {
165  const double y_dbl = value_of(y_vec[i]);
166  const double mu_dbl = value_of(mu_vec[i]);
167  const double sigma_dbl = value_of(sigma_vec[i]);
168  const double nu_dbl = value_of(nu_vec[i]);
169  square_y_minus_mu_over_sigma__over_nu[i]
170  = square((y_dbl - mu_dbl) / sigma_dbl) / nu_dbl;
171  log1p_exp[i] = log1p(square_y_minus_mu_over_sigma__over_nu[i]);
172  }
173 
174  agrad::OperandsAndPartials<T_y,T_dof,T_loc,T_scale>
175  operands_and_partials(y,nu,mu,sigma);
176  for (size_t n = 0; n < N; n++) {
177  const double y_dbl = value_of(y_vec[n]);
178  const double mu_dbl = value_of(mu_vec[n]);
179  const double sigma_dbl = value_of(sigma_vec[n]);
180  const double nu_dbl = value_of(nu_vec[n]);
181  if (include_summand<propto>::value)
182  logp += NEG_LOG_SQRT_PI;
183  if (include_summand<propto,T_dof>::value)
184  logp += lgamma_half_nu_plus_half[n] - lgamma_half_nu[n]
185  - 0.5 * log_nu[n];
186  if (include_summand<propto,T_scale>::value)
187  logp -= log_sigma[n];
188  if (include_summand<propto,T_y,T_dof,T_loc,T_scale>::value)
189  logp -= (half_nu[n] + 0.5)
190  * log1p_exp[n];
191 
192  if (!is_constant_struct<T_y>::value) {
193  operands_and_partials.d_x1[n]
194  += -(half_nu[n]+0.5)
195  * 1.0 / (1.0 + square_y_minus_mu_over_sigma__over_nu[n])
196  * (2.0 * (y_dbl - mu_dbl) / square(sigma_dbl) / nu_dbl);
197  }
198  if (!is_constant_struct<T_dof>::value) {
199  const double inv_nu = 1.0 / nu_dbl;
200  operands_and_partials.d_x2[n]
201  += 0.5*digamma_half_nu_plus_half[n] - 0.5*digamma_half_nu[n]
202  - 0.5 * inv_nu
203  - 0.5*log1p_exp[n]
204  + (half_nu[n] + 0.5)
205  * (1.0/(1.0 + square_y_minus_mu_over_sigma__over_nu[n])
206  * square_y_minus_mu_over_sigma__over_nu[n] * inv_nu);
207  }
208  if (!is_constant_struct<T_loc>::value) {
209  operands_and_partials.d_x3[n]
210  -= (half_nu[n] + 0.5)
211  / (1.0 + square_y_minus_mu_over_sigma__over_nu[n])
212  * (2.0 * (mu_dbl - y_dbl) / (sigma_dbl*sigma_dbl*nu_dbl));
213  }
214  if (!is_constant_struct<T_scale>::value) {
215  const double inv_sigma = 1.0 / sigma_dbl;
216  operands_and_partials.d_x4[n]
217  += -inv_sigma
218  + (nu_dbl + 1.0) / (1.0 + square_y_minus_mu_over_sigma__over_nu[n])
219  * (square_y_minus_mu_over_sigma__over_nu[n] * inv_sigma);
220  }
221  }
222  return operands_and_partials.to_var(logp);
223  }
224 
225  template <bool propto,
226  typename T_y, typename T_dof, typename T_loc, typename T_scale>
227  inline
228  typename return_type<T_y,T_dof,T_loc,T_scale>::type
229  student_t_log(const T_y& y, const T_dof& nu, const T_loc& mu,
230  const T_scale& sigma) {
231  return student_t_log<propto>(y,nu,mu,sigma,stan::math::default_policy());
232  }
233 
234  template <typename T_y, typename T_dof, typename T_loc, typename T_scale,
235  class Policy>
236  inline
237  typename return_type<T_y,T_dof,T_loc,T_scale>::type
238  student_t_log(const T_y& y, const T_dof& nu, const T_loc& mu,
239  const T_scale& sigma,
240  const Policy&) {
241  return student_t_log<false>(y,nu,mu,sigma,Policy());
242  }
243 
244  template <typename T_y, typename T_dof, typename T_loc, typename T_scale>
245  inline
246  typename return_type<T_y,T_dof,T_loc,T_scale>::type
247  student_t_log(const T_y& y, const T_dof& nu, const T_loc& mu,
248  const T_scale& sigma) {
249  return student_t_log<false>(y,nu,mu,sigma,stan::math::default_policy());
250  }
251 
252  template <typename T_y, typename T_dof, typename T_loc, typename T_scale,
253  class Policy>
254  typename return_type<T_y, T_dof, T_loc, T_scale>::type
255  student_t_cdf(const T_y& y, const T_dof& nu, const T_loc& mu,
256  const T_scale& sigma, const Policy&) {
257 
258  // Size checks
259  if (!(stan::length(y) && stan::length(nu) && stan::length(mu)
260  && stan::length(sigma)))
261  return 1.0;
262 
263  static const char* function = "stan::prob::student_t_cdf(%1%)";
264 
265  using stan::math::check_positive;
266  using stan::math::check_finite;
267  using stan::math::check_not_nan;
268  using stan::math::check_consistent_sizes;
269 
270  using stan::math::value_of;
271 
272  double P(1.0);
273 
274  if (!check_not_nan(function, y, "Random variable", &P, Policy()))
275  return P;
276 
277  if (!check_finite(function, nu, "Degrees of freedom parameter", &P,
278  Policy()))
279  return P;
280 
281  if (!check_positive(function, nu, "Degrees of freedom parameter", &P,
282  Policy()))
283  return P;
284 
285  if (!check_finite(function, mu, "Location parameter", &P, Policy()))
286  return P;
287 
288  if (!check_finite(function, sigma, "Scale parameter", &P, Policy()))
289  return P;
290 
291  if (!check_positive(function, sigma, "Scale parameter", &P, Policy()))
292  return P;
293 
294  // Wrap arguments in vectors
295  VectorView<const T_y> y_vec(y);
296  VectorView<const T_dof> nu_vec(nu);
297  VectorView<const T_loc> mu_vec(mu);
298  VectorView<const T_scale> sigma_vec(sigma);
299  size_t N = max_size(y, nu, mu, sigma);
300 
301  agrad::OperandsAndPartials<T_y, T_dof, T_loc, T_scale>
302  operands_and_partials(y, nu, mu, sigma);
303 
304  std::fill(operands_and_partials.all_partials,
305  operands_and_partials.all_partials
306  + operands_and_partials.nvaris,
307  0.0);
308 
309  // Explicit return for extreme values
310  // The gradients are technically ill-defined, but treated as zero
311  for (size_t i = 0; i < stan::length(y); i++) {
312  if (value_of(y_vec[i]) == -std::numeric_limits<double>::infinity())
313  return operands_and_partials.to_var(0.0);
314  }
315 
316  using boost::math::ibeta;
317  using boost::math::ibeta_derivative;
318 
319  using boost::math::digamma;
320  using boost::math::beta;
321 
322  // Cache a few expensive function calls if nu is a parameter
323  double digammaHalf = 0;
324 
325  DoubleVectorView<!is_constant_struct<T_dof>::value,is_vector<T_dof>::value>
326  digamma_vec(stan::length(nu));
327  DoubleVectorView<!is_constant_struct<T_dof>::value,is_vector<T_dof>::value>
328  digammaNu_vec(stan::length(nu));
329  DoubleVectorView<!is_constant_struct<T_dof>::value,is_vector<T_dof>::value>
330  digammaNuPlusHalf_vec(stan::length(nu));
331  DoubleVectorView<!is_constant_struct<T_dof>::value,is_vector<T_dof>::value>
332  betaNuHalf_vec(stan::length(nu));
333 
334  if (!is_constant_struct<T_dof>::value) {
335 
336  digammaHalf = digamma(0.5);
337 
338  for (size_t i = 0; i < stan::length(nu); i++) {
339 
340  const double nu_dbl = value_of(nu_vec[i]);
341 
342  digammaNu_vec[i] = digamma(0.5 * nu_dbl);
343  digammaNuPlusHalf_vec[i] = digamma(0.5 + 0.5 * nu_dbl);
344  betaNuHalf_vec[i] = beta(0.5, 0.5 * nu_dbl);
345 
346  }
347 
348  }
349 
350  // Compute vectorized CDF and gradient
351  for (size_t n = 0; n < N; n++) {
352 
353  // Explicit results for extreme values
354  // The gradients are technically ill-defined, but treated as zero
355  if (value_of(y_vec[n]) == std::numeric_limits<double>::infinity()) {
356  continue;
357  }
358 
359  const double sigma_inv = 1.0 / value_of(sigma_vec[n]);
360  const double t = (value_of(y_vec[n]) - value_of(mu_vec[n])) * sigma_inv;
361  const double nu_dbl = value_of(nu_vec[n]);
362  const double q = nu_dbl / (t * t);
363  const double r = 1.0 / (1.0 + q);
364  const double J = 2 * r * r * q / t;
365  double zJacobian = t > 0 ? - 0.5 : 0.5;
366 
367  if(q < 2)
368  {
369 
370  double z = ibeta(0.5 * nu_dbl, 0.5, 1.0 - r);
371 
372  const double Pn = t > 0 ? 1.0 - 0.5 * z : 0.5 * z;
373 
374  const double d_ibeta = ibeta_derivative(0.5 * nu_dbl, 0.5, 1.0 - r);
375 
376  P *= Pn;
377 
378  if (!is_constant_struct<T_y>::value)
379  operands_and_partials.d_x1[n]
380  += - zJacobian * d_ibeta * J * sigma_inv / Pn;
381 
382  if (!is_constant_struct<T_dof>::value) {
383 
384  double g1 = 0;
385  double g2 = 0;
386 
387  stan::math::gradRegIncBeta(g1, g2, 0.5 * nu_dbl, 0.5, 1.0 - r,
388  digammaNu_vec[n], digammaHalf,
389  digammaNuPlusHalf_vec[n],
390  betaNuHalf_vec[n]);
391 
392  operands_and_partials.d_x2[n]
393  += zJacobian * ( d_ibeta * (r / t) * (r / t) + 0.5 * g1 ) / Pn;
394 
395  }
396 
397  if (!is_constant_struct<T_loc>::value)
398  operands_and_partials.d_x3[n]
399  += zJacobian * d_ibeta * J * sigma_inv / Pn;
400 
401  if (!is_constant_struct<T_scale>::value)
402  operands_and_partials.d_x4[n]
403  += zJacobian * d_ibeta * J * sigma_inv * t / Pn;
404 
405  }
406  else
407  {
408 
409  double z = 1 - ibeta(0.5, 0.5 * nu_dbl, r);
410  zJacobian *= -1;
411 
412  const double Pn = t > 0 ? 1.0 - 0.5 * z : 0.5 * z;
413 
414  double d_ibeta = ibeta_derivative(0.5, 0.5 * nu_dbl, r);
415 
416  P *= Pn;
417 
418  if (!is_constant_struct<T_y>::value)
419  operands_and_partials.d_x1[n]
420  += zJacobian * d_ibeta * J * sigma_inv / Pn;
421 
422  if (!is_constant_struct<T_dof>::value) {
423 
424  double g1 = 0;
425  double g2 = 0;
426 
427  stan::math::gradRegIncBeta(g1, g2, 0.5, 0.5 * nu_dbl, r,
428  digammaHalf, digammaNu_vec[n],
429  digammaNuPlusHalf_vec[n],
430  betaNuHalf_vec[n]);
431 
432  operands_and_partials.d_x2[n]
433  += zJacobian * ( - d_ibeta * (r / t) * (r / t) + 0.5 * g2 ) / Pn;
434 
435  }
436 
437  if (!is_constant_struct<T_loc>::value)
438  operands_and_partials.d_x3[n]
439  += - zJacobian * d_ibeta * J * sigma_inv / Pn;
440 
441  if (!is_constant_struct<T_scale>::value)
442  operands_and_partials.d_x4[n]
443  += - zJacobian * d_ibeta * J * sigma_inv * t / Pn;
444  }
445 
446  }
447 
448  if (!is_constant_struct<T_y>::value)
449  for(size_t n = 0; n < stan::length(y); ++n)
450  operands_and_partials.d_x1[n] *= P;
451 
452  if (!is_constant_struct<T_dof>::value)
453  for(size_t n = 0; n < stan::length(nu); ++n)
454  operands_and_partials.d_x2[n] *= P;
455 
456  if (!is_constant_struct<T_loc>::value)
457  for(size_t n = 0; n < stan::length(mu); ++n)
458  operands_and_partials.d_x3[n] *= P;
459 
460 
461  if (!is_constant_struct<T_scale>::value)
462  for(size_t n = 0; n < stan::length(sigma); ++n)
463  operands_and_partials.d_x4[n] *= P;
464 
465 
466  return operands_and_partials.to_var(P);
467 
468  }
469 
470  template <typename T_y, typename T_dof, typename T_loc, typename T_scale>
471  typename return_type<T_y,T_dof,T_loc,T_scale>::type
472  student_t_cdf(const T_y& y, const T_dof& nu, const T_loc& mu, const T_scale& sigma) {
473  return student_t_cdf(y, nu, mu, sigma, stan::math::default_policy());
474  }
475 
476 
477  template <class RNG>
478  inline double
479  student_t_rng(const double nu,
480  const double mu,
481  const double sigma,
482  RNG& rng) {
483  using boost::variate_generator;
484  using boost::random::student_t_distribution;
485  variate_generator<RNG&, student_t_distribution<> >
486  rng_unit_student_t(rng, student_t_distribution<>(nu));
487  return mu + sigma * rng_unit_student_t();
488  }
489  }
490 }
491 #endif
     [ Stan Home Page ] © 2011–2013, Stan Development Team.