Stan  1.1
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 <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 
35  template <bool propto,
36  typename T_y, typename T_dof, typename T_scale,
37  class Policy>
38  typename return_type<T_y,T_dof,T_scale>::type
39  scaled_inv_chi_square_log(const T_y& y, const T_dof& nu, const T_scale& s,
40  const Policy&) {
41  static const char* function
42  = "stan::prob::scaled_inv_chi_square_log(%1%)";
43 
44  using stan::math::check_finite;
45  using stan::math::check_positive;
46  using stan::math::check_not_nan;
47  using stan::math::check_consistent_sizes;
48  using stan::math::value_of;
49 
50  // check if any vectors are zero length
51  if (!(stan::length(y)
52  && stan::length(nu)
53  && stan::length(s)))
54  return 0.0;
55 
56  typename return_type<T_y,T_dof,T_scale>::type logp(0.0);
57  if (!check_not_nan(function, y, "Random variable", &logp, Policy()))
58  return logp;
59  if (!check_finite(function, nu, "Degrees of freedom parameter", &logp, Policy()))
60  return logp;
61  if (!check_positive(function, nu, "Degrees of freedom parameter", &logp, Policy()))
62  return logp;
63  if (!check_finite(function, s, "Scale parameter", &logp, Policy()))
64  return logp;
65  if (!check_positive(function, s, "Scale parameter", &logp, Policy()))
66  return logp;
67  if (!(check_consistent_sizes(function,
68  y,nu,s,
69  "Random variable","Degrees of freedom parameter","Scale parameter",
70  &logp, Policy())))
71  return logp;
72 
73  // check if no variables are involved and prop-to
74  if (!include_summand<propto,T_y,T_dof,T_scale>::value)
75  return 0.0;
76 
77  VectorView<const T_y> y_vec(y);
78  VectorView<const T_dof> nu_vec(nu);
79  VectorView<const T_scale> s_vec(s);
80  size_t N = max_size(y, nu, s);
81 
82  for (size_t n = 0; n < N; n++) {
83  if (value_of(y_vec[n]) <= 0)
84  return LOG_ZERO;
85  }
86 
87  using stan::math::multiply_log;
88  using stan::math::square;
89 
90  for (size_t n = 0; n < N; n++) {
91  if (include_summand<propto,T_dof>::value) {
92  typename return_type<T_dof>::type half_nu = 0.5 * nu_vec[n];
93  logp += multiply_log(half_nu,half_nu) - lgamma(half_nu);
94  }
95  if (include_summand<propto,T_dof,T_scale>::value)
96  logp += nu_vec[n] * log(s_vec[n]);
97  if (include_summand<propto,T_dof,T_y>::value)
98  logp -= multiply_log(nu_vec[n]*0.5+1.0, y_vec[n]);
99  if (include_summand<propto,T_dof,T_y,T_scale>::value)
100  logp -= nu_vec[n] * 0.5 * square(s_vec[n]) / y_vec[n];
101  }
102  return logp;
103  }
104 
105  template <bool propto,
106  typename T_y, typename T_dof, typename T_scale>
107  inline
108  typename return_type<T_y,T_dof,T_scale>::type
109  scaled_inv_chi_square_log(const T_y& y, const T_dof& nu, const T_scale& s) {
110  return scaled_inv_chi_square_log<propto>(y,nu,s,
111  stan::math::default_policy());
112  }
113 
114  template <typename T_y, typename T_dof, typename T_scale,
115  class Policy>
116  inline
117  typename return_type<T_y,T_dof,T_scale>::type
118  scaled_inv_chi_square_log(const T_y& y, const T_dof& nu, const T_scale& s,
119  const Policy&) {
120  return scaled_inv_chi_square_log<false>(y,nu,s,Policy());
121  }
122 
123  template <typename T_y, typename T_dof, typename T_scale>
124  inline
125  typename return_type<T_y,T_dof,T_scale>::type
126  scaled_inv_chi_square_log(const T_y& y, const T_dof& nu, const T_scale& s) {
127  return scaled_inv_chi_square_log<false>(y,nu,s,
128  stan::math::default_policy());
129  }
130 
145  template <typename T_y, typename T_dof, typename T_scale, class Policy>
146  typename return_type<T_y, T_dof, T_scale>::type
147  scaled_inv_chi_square_cdf(const T_y& y, const T_dof& nu, const T_scale& s, const Policy&) {
148 
149  static const char* function = "stan::prob::scaled_inv_chi_square_log(%1%)";
150 
151  using stan::math::check_finite;
152  using stan::math::check_positive;
153  using stan::math::check_not_nan;
154  using stan::math::check_consistent_sizes;
155 
156  using boost::math::tools::promote_args;
157 
158  double P(1.0);
159 
160  if (!check_not_nan(function, y, "Random variable", &P, Policy()))
161  return P;
162 
163  if (!check_positive(function, y, "Random variable", &P, Policy()))
164  return P;
165 
166  if (!check_finite(function, nu, "Degrees of freedom parameter", &P, Policy()))
167  return P;
168 
169  if (!check_positive(function, nu, "Degrees of freedom parameter", &P, Policy()))
170  return P;
171 
172  if (!check_finite(function, s, "Scale parameter", &P, Policy()))
173  return P;
174 
175  if (!check_positive(function, s, "Scale parameter", &P, Policy()))
176  return P;
177 
178  if (!(check_consistent_sizes(function, y, nu, s,
179  "Random variable", "Degrees of freedom 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_dof> nu_vec(nu);
186  VectorView<const T_scale> s_vec(s);
187  size_t N = max_size(y, nu, s);
188 
189  agrad::OperandsAndPartials<T_y, T_dof, T_scale> operands_and_partials(y, nu, s);
190 
191  std::fill(operands_and_partials.all_partials,
192  operands_and_partials.all_partials + operands_and_partials.nvaris, 0.0);
193 
194  // Compute CDF and its gradients
195  using stan::math::value_of;
196  using boost::math::gamma_p_derivative;
197  using boost::math::gamma_q;
198  using boost::math::digamma;
199  using boost::math::tgamma;
200 
201  // Cache a few expensive function calls if nu is a parameter
202  DoubleVectorView<!is_constant_struct<T_dof>::value,is_vector<T_dof>::value> gamma_vec(stan::length(nu));
203  DoubleVectorView<!is_constant_struct<T_dof>::value,is_vector<T_dof>::value> digamma_vec(stan::length(nu));
204 
205  if (!is_constant_struct<T_dof>::value) {
206 
207  for (size_t i = 0; i < stan::length(nu); i++) {
208  const double half_nu_dbl = 0.5 * value_of(nu_vec[i]);
209  gamma_vec[i] = tgamma(half_nu_dbl);
210  digamma_vec[i] = digamma(half_nu_dbl);
211  }
212 
213  }
214 
215  // Compute vectorized CDF and gradient
216  for (size_t n = 0; n < N; n++) {
217 
218  // Pull out values
219  const double y_dbl = value_of(y_vec[n]);
220  const double y_inv_dbl = 1.0 / y_dbl;
221  const double half_nu_dbl = 0.5 * value_of(nu_vec[n]);
222  const double s_dbl = value_of(s_vec[n]);
223  const double half_s2_overx_dbl = 0.5 * s_dbl * s_dbl * y_inv_dbl;
224  const double half_nu_s2_overx_dbl = 2.0 * half_nu_dbl * half_s2_overx_dbl;
225 
226  // Compute
227  const double Pn = gamma_q(half_nu_dbl, half_nu_s2_overx_dbl);
228 
229  P *= Pn;
230 
231  if (!is_constant_struct<T_y>::value)
232  operands_and_partials.d_x1[n]
233  += half_nu_s2_overx_dbl * y_inv_dbl * gamma_p_derivative(half_nu_dbl, half_nu_s2_overx_dbl) / Pn;
234 
235  if (!is_constant_struct<T_dof>::value)
236  operands_and_partials.d_x2[n]
237  += (0.5 * stan::math::gradRegIncGamma(half_nu_dbl, half_nu_s2_overx_dbl, gamma_vec[n], digamma_vec[n])
238  - half_s2_overx_dbl * gamma_p_derivative(half_nu_dbl, half_nu_s2_overx_dbl) ) / Pn;
239 
240  if (!is_constant_struct<T_scale>::value)
241  operands_and_partials.d_x3[n]
242  += - 2.0 * half_nu_dbl * s_dbl * y_inv_dbl * gamma_p_derivative(half_nu_dbl, half_nu_s2_overx_dbl) / Pn;
243 
244  }
245 
246  for (size_t n = 0; n < N; n++) {
247 
248  if (!is_constant_struct<T_y>::value)
249  operands_and_partials.d_x1[n] *= P;
250 
251  if (!is_constant_struct<T_dof>::value)
252  operands_and_partials.d_x2[n] *= P;
253 
254  if (!is_constant_struct<T_scale>::value)
255  operands_and_partials.d_x3[n] *= P;
256 
257  }
258 
259  return operands_and_partials.to_var(P);
260 
261  }
262 
263 
264  template <typename T_y, typename T_dof, typename T_scale>
265  inline typename return_type<T_y, T_dof, T_scale>::type
266  scaled_inv_chi_square_cdf(const T_y& y, const T_dof& nu, const T_scale& s) {
267  return scaled_inv_chi_square_cdf(y, nu, s, stan::math::default_policy());
268  }
269 
270  }
271 }
272 
273 #endif
274 
     [ Stan Home Page ] © 2011–2013, Stan Development Team.