Stan  1.3
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
scaled_inv_chi_square.hpp
Go to the documentation of this file.
1 #ifndef __STAN__PROB__DISTRIBUTIONS__UNIVARIATE__CONTINUOUS__SCALED_INV_CHI_SQUARE_HPP__
2 #define __STAN__PROB__DISTRIBUTIONS__UNIVARIATE__CONTINUOUS__SCALED_INV_CHI_SQUARE_HPP__
3 
4 #include <boost/random/chi_squared_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 
39  template <bool propto,
40  typename T_y, typename T_dof, typename T_scale,
41  class Policy>
42  typename return_type<T_y,T_dof,T_scale>::type
43  scaled_inv_chi_square_log(const T_y& y, const T_dof& nu, const T_scale& s,
44  const Policy&) {
45  static const char* function
46  = "stan::prob::scaled_inv_chi_square_log(%1%)";
47 
48  using stan::math::check_finite;
49  using stan::math::check_positive;
50  using stan::math::check_not_nan;
51  using stan::math::check_consistent_sizes;
52  using stan::math::value_of;
53 
54  // check if any vectors are zero length
55  if (!(stan::length(y)
56  && stan::length(nu)
57  && stan::length(s)))
58  return 0.0;
59 
60  double logp(0.0);
61  if (!check_not_nan(function, y, "Random variable", &logp, Policy()))
62  return logp;
63  if (!check_finite(function, nu, "Degrees of freedom parameter",
64  &logp, Policy()))
65  return logp;
66  if (!check_positive(function, nu, "Degrees of freedom parameter",
67  &logp, Policy()))
68  return logp;
69  if (!check_finite(function, s, "Scale parameter", &logp, Policy()))
70  return logp;
71  if (!check_positive(function, s, "Scale parameter", &logp, Policy()))
72  return logp;
73  if (!(check_consistent_sizes(function,
74  y,nu,s,
75  "Random variable",
76  "Degrees of freedom parameter",
77  "Scale parameter",
78  &logp, Policy())))
79  return logp;
80 
81  // check if no variables are involved and prop-to
82  if (!include_summand<propto,T_y,T_dof,T_scale>::value)
83  return 0.0;
84 
85  VectorView<const T_y> y_vec(y);
86  VectorView<const T_dof> nu_vec(nu);
87  VectorView<const T_scale> s_vec(s);
88  size_t N = max_size(y, nu, s);
89 
90  for (size_t n = 0; n < N; n++) {
91  if (value_of(y_vec[n]) <= 0)
92  return LOG_ZERO;
93  }
94 
95  using boost::math::lgamma;
96  using boost::math::digamma;
97  using stan::math::multiply_log;
98  using stan::math::square;
99 
100  DoubleVectorView<include_summand<propto,T_dof,T_y,T_scale>::value,
101  is_vector<T_dof>::value> half_nu(length(nu));
102  for (size_t i = 0; i < length(nu); i++)
103  if (include_summand<propto,T_dof,T_y,T_scale>::value)
104  half_nu[i] = 0.5 * value_of(nu_vec[i]);
105 
106  DoubleVectorView<include_summand<propto,T_dof,T_y>::value,
107  is_vector<T_y>::value> log_y(length(y));
108  for (size_t i = 0; i < length(y); i++)
109  if (include_summand<propto,T_dof,T_y>::value)
110  log_y[i] = log(value_of(y_vec[i]));
111 
112  DoubleVectorView<include_summand<propto,T_dof,T_y,T_scale>::value,
113  is_vector<T_y>::value> inv_y(length(y));
114  for (size_t i = 0; i < length(y); i++)
115  if (include_summand<propto,T_dof,T_y,T_scale>::value)
116  inv_y[i] = 1.0 / value_of(y_vec[i]);
117 
118  DoubleVectorView<include_summand<propto,T_dof,T_scale>::value,
119  is_vector<T_scale>::value> log_s(length(s));
120  for (size_t i = 0; i < length(s); i++)
121  if (include_summand<propto,T_dof,T_scale>::value)
122  log_s[i] = log(value_of(s_vec[i]));
123 
124  DoubleVectorView<include_summand<propto,T_dof>::value,
125  is_vector<T_dof>::value> log_half_nu(length(nu));
126  DoubleVectorView<include_summand<propto,T_dof>::value,
127  is_vector<T_dof>::value> lgamma_half_nu(length(nu));
128  DoubleVectorView<!is_constant_struct<T_dof>::value,
129  is_vector<T_dof>::value> digamma_half_nu_over_two(length(nu));
130  for (size_t i = 0; i < length(nu); i++) {
131  if (include_summand<propto,T_dof>::value)
132  lgamma_half_nu[i] = lgamma(half_nu[i]);
133  if (include_summand<propto,T_dof>::value)
134  log_half_nu[i] = log(half_nu[i]);
135  if (!is_constant_struct<T_dof>::value)
136  digamma_half_nu_over_two[i] = digamma(half_nu[i]) * 0.5;
137  }
138 
139  agrad::OperandsAndPartials<T_y,T_dof,T_scale>
140  operands_and_partials(y, nu, s);
141  for (size_t n = 0; n < N; n++) {
142  const double s_dbl = value_of(s_vec[n]);
143  const double nu_dbl = value_of(nu_vec[n]);
144  if (include_summand<propto,T_dof>::value)
145  logp += half_nu[n] * log_half_nu[n] - lgamma_half_nu[n];
146  if (include_summand<propto,T_dof,T_scale>::value)
147  logp += nu_dbl * log_s[n];
148  if (include_summand<propto,T_dof,T_y>::value)
149  logp -= (half_nu[n]+1.0) * log_y[n];
150  if (include_summand<propto,T_dof,T_y,T_scale>::value)
151  logp -= half_nu[n] * s_dbl*s_dbl * inv_y[n];
152 
153  if (!is_constant_struct<T_y>::value) {
154  operands_and_partials.d_x1[n]
155  += -(half_nu[n] + 1.0) * inv_y[n]
156  + half_nu[n] * s_dbl*s_dbl * inv_y[n]*inv_y[n];
157  }
158  if (!is_constant_struct<T_dof>::value) {
159  operands_and_partials.d_x2[n]
160  += 0.5 * log_half_nu[n] + 0.5
161  - digamma_half_nu_over_two[n]
162  + log_s[n]
163  - 0.5 * log_y[n]
164  - 0.5* s_dbl*s_dbl * inv_y[n];
165  }
166  if (!is_constant_struct<T_scale>::value) {
167  operands_and_partials.d_x3[n]
168  += nu_dbl / s_dbl - nu_dbl * inv_y[n] * s_dbl;
169  }
170  }
171  return operands_and_partials.to_var(logp);
172  }
173 
174  template <bool propto,
175  typename T_y, typename T_dof, typename T_scale>
176  inline
177  typename return_type<T_y,T_dof,T_scale>::type
178  scaled_inv_chi_square_log(const T_y& y, const T_dof& nu,
179  const T_scale& s) {
180  return scaled_inv_chi_square_log<propto>(y,nu,s,
181  stan::math::default_policy());
182  }
183 
184  template <typename T_y, typename T_dof, typename T_scale,
185  class Policy>
186  inline
187  typename return_type<T_y,T_dof,T_scale>::type
188  scaled_inv_chi_square_log(const T_y& y, const T_dof& nu, const T_scale& s,
189  const Policy&) {
190  return scaled_inv_chi_square_log<false>(y,nu,s,Policy());
191  }
192 
193  template <typename T_y, typename T_dof, typename T_scale>
194  inline
195  typename return_type<T_y,T_dof,T_scale>::type
196  scaled_inv_chi_square_log(const T_y& y, const T_dof& nu, const T_scale& s) {
197  return scaled_inv_chi_square_log<false>(y,nu,s,
198  stan::math::default_policy());
199  }
200 
215  template <typename T_y, typename T_dof, typename T_scale, class Policy>
216  typename return_type<T_y, T_dof, T_scale>::type
217  scaled_inv_chi_square_cdf(const T_y& y, const T_dof& nu,
218  const T_scale& s, const Policy&) {
219  // Size checks
220  if (!(stan::length(y) && stan::length(nu) && stan::length(s)))
221  return 1.0;
222 
223  static const char* function
224  = "stan::prob::scaled_inv_chi_square_log(%1%)";
225 
226  using stan::math::check_finite;
227  using stan::math::check_positive;
228  using stan::math::check_not_nan;
229  using stan::math::check_consistent_sizes;
230  using stan::math::check_nonnegative;
231  using stan::math::value_of;
232 
233  double P(1.0);
234 
235  if (!check_not_nan(function, y, "Random variable", &P, Policy()))
236  return P;
237 
238  if (!check_nonnegative(function, y, "Random variable", &P, Policy()))
239  return P;
240 
241  if (!check_finite(function, nu, "Degrees of freedom parameter",
242  &P, Policy()))
243  return P;
244 
245  if (!check_positive(function, nu, "Degrees of freedom parameter",
246  &P, Policy()))
247  return P;
248 
249  if (!check_finite(function, s, "Scale parameter", &P, Policy()))
250  return P;
251 
252  if (!check_positive(function, s, "Scale parameter", &P, Policy()))
253  return P;
254 
255  if (!(check_consistent_sizes(function, y, nu, s,
256  "Random variable",
257  "Degrees of freedom parameter",
258  "Scale parameter",
259  &P, Policy())))
260  return P;
261 
262  // Wrap arguments in vectors
263  VectorView<const T_y> y_vec(y);
264  VectorView<const T_dof> nu_vec(nu);
265  VectorView<const T_scale> s_vec(s);
266  size_t N = max_size(y, nu, s);
267 
268  agrad::OperandsAndPartials<T_y, T_dof, T_scale>
269  operands_and_partials(y, nu, s);
270 
271  std::fill(operands_and_partials.all_partials,
272  operands_and_partials.all_partials
273  + operands_and_partials.nvaris, 0.0);
274 
275  // Explicit return for extreme values
276  // The gradients are technically ill-defined, but treated as zero
277 
278  for (size_t i = 0; i < stan::length(y); i++) {
279  if (value_of(y_vec[i]) == 0)
280  return operands_and_partials.to_var(0.0);
281  }
282 
283  // Compute CDF and its gradients
284  using boost::math::gamma_p_derivative;
285  using boost::math::gamma_q;
286  using boost::math::digamma;
287  using boost::math::tgamma;
288 
289  // Cache a few expensive function calls if nu is a parameter
290  DoubleVectorView<!is_constant_struct<T_dof>::value,
291  is_vector<T_dof>::value> gamma_vec(stan::length(nu));
292  DoubleVectorView<!is_constant_struct<T_dof>::value,
293  is_vector<T_dof>::value> digamma_vec(stan::length(nu));
294 
295  if (!is_constant_struct<T_dof>::value) {
296 
297  for (size_t i = 0; i < stan::length(nu); i++) {
298  const double half_nu_dbl = 0.5 * value_of(nu_vec[i]);
299  gamma_vec[i] = tgamma(half_nu_dbl);
300  digamma_vec[i] = digamma(half_nu_dbl);
301  }
302 
303  }
304 
305  // Compute vectorized CDF and gradient
306  for (size_t n = 0; n < N; n++) {
307 
308  // Explicit results for extreme values
309  // The gradients are technically ill-defined, but treated as zero
310  if (value_of(y_vec[n]) == std::numeric_limits<double>::infinity()) {
311  continue;
312  }
313 
314  // Pull out values
315  const double y_dbl = value_of(y_vec[n]);
316  const double y_inv_dbl = 1.0 / y_dbl;
317  const double half_nu_dbl = 0.5 * value_of(nu_vec[n]);
318  const double s_dbl = value_of(s_vec[n]);
319  const double half_s2_overx_dbl = 0.5 * s_dbl * s_dbl * y_inv_dbl;
320  const double half_nu_s2_overx_dbl
321  = 2.0 * half_nu_dbl * half_s2_overx_dbl;
322 
323  // Compute
324  const double Pn = gamma_q(half_nu_dbl, half_nu_s2_overx_dbl);
325 
326  P *= Pn;
327 
328  if (!is_constant_struct<T_y>::value)
329  operands_and_partials.d_x1[n]
330  += half_nu_s2_overx_dbl * y_inv_dbl
331  * gamma_p_derivative(half_nu_dbl, half_nu_s2_overx_dbl) / Pn;
332 
333  if (!is_constant_struct<T_dof>::value)
334  operands_and_partials.d_x2[n]
335  += (0.5 * stan::math::gradRegIncGamma(half_nu_dbl,
336  half_nu_s2_overx_dbl,
337  gamma_vec[n], digamma_vec[n])
338  - half_s2_overx_dbl
339  * gamma_p_derivative(half_nu_dbl, half_nu_s2_overx_dbl) )
340  / Pn;
341 
342  if (!is_constant_struct<T_scale>::value)
343  operands_and_partials.d_x3[n]
344  += - 2.0 * half_nu_dbl * s_dbl * y_inv_dbl
345  * gamma_p_derivative(half_nu_dbl, half_nu_s2_overx_dbl) / Pn;
346 
347  }
348 
349  if (!is_constant_struct<T_y>::value)
350  for(size_t n = 0; n < stan::length(y); ++n)
351  operands_and_partials.d_x1[n] *= P;
352 
353  if (!is_constant_struct<T_dof>::value)
354  for(size_t n = 0; n < stan::length(nu); ++n)
355  operands_and_partials.d_x2[n] *= P;
356 
357  if (!is_constant_struct<T_scale>::value)
358  for(size_t n = 0; n < stan::length(s); ++n)
359  operands_and_partials.d_x3[n] *= P;
360 
361  return operands_and_partials.to_var(P);
362 
363  }
364 
365 
366  template <typename T_y, typename T_dof, typename T_scale>
367  inline typename return_type<T_y, T_dof, T_scale>::type
368  scaled_inv_chi_square_cdf(const T_y& y, const T_dof& nu,
369  const T_scale& s) {
370  return scaled_inv_chi_square_cdf(y, nu, s,
371  stan::math::default_policy());
372  }
373 
374  template <class RNG>
375  inline double
376  scaled_inv_chi_square_rng(const double nu,
377  const double s,
378  RNG& rng) {
379  using boost::variate_generator;
380  using boost::random::chi_squared_distribution;
381  variate_generator<RNG&, chi_squared_distribution<> >
382  chi_square_rng(rng, chi_squared_distribution<>(nu));
383  return nu * s / chi_square_rng();
384  }
385  }
386 }
387 
388 #endif
389 
     [ Stan Home Page ] © 2011–2013, Stan Development Team.