Stan  1.0
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
inv_gamma.hpp
Go to the documentation of this file.
1 #ifndef __STAN__PROB__DISTRIBUTIONS__UNIVARIATE__CONTINUOUS__INV_GAMMA_HPP__
2 #define __STAN__PROB__DISTRIBUTIONS__UNIVARIATE__CONTINUOUS__INV_GAMMA_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 
11 namespace stan {
12 
13  namespace prob {
14 
31  template <bool propto,
32  typename T_y, typename T_shape, typename T_scale,
33  class Policy>
34  typename return_type<T_y,T_shape,T_scale>::type
35  inv_gamma_log(const T_y& y, const T_shape& alpha, const T_scale& beta,
36  const Policy&) {
37  static const char* function = "stan::prob::inv_gamma_log(%1%)";
38 
39  using stan::is_constant_struct;
40  using stan::math::check_not_nan;
41  using stan::math::check_positive;
42  using stan::math::check_finite;
43  using boost::math::tools::promote_args;
44  using stan::math::check_consistent_sizes;
45  using stan::math::value_of;
46 
47  // check if any vectors are zero length
48  if (!(stan::length(y)
49  && stan::length(alpha)
50  && stan::length(beta)))
51  return 0.0;
52 
53  // set up return value accumulator
54  double logp(0.0);
55 
56  if (!check_not_nan(function, y, "Random variable", &logp, Policy()))
57  return logp;
58  if (!check_finite(function, alpha, "Shape parameter",
59  &logp, Policy()))
60  return logp;
61  if (!check_positive(function, alpha, "Shape parameter",
62  &logp, Policy()))
63  return logp;
64  if (!check_finite(function, beta, "Scale parameter",
65  &logp, Policy()))
66  return logp;
67  if (!check_positive(function, beta, "Scale parameter",
68  &logp, Policy()))
69  return logp;
70  if (!(check_consistent_sizes(function,
71  y,alpha,beta,
72  "Random variable","Shape parameter","Scale parameter",
73  &logp, Policy())))
74  return logp;
75 
76  // check if no variables are involved and prop-to
77  if (!include_summand<propto,T_y,T_shape,T_scale>::value)
78  return 0.0;
79 
80  // set up template expressions wrapping scalars into vector views
81  VectorView<const T_y> y_vec(y);
82  VectorView<const T_shape> alpha_vec(alpha);
83  VectorView<const T_scale> beta_vec(beta);
84 
85  for (size_t n = 0; n < length(y); n++) {
86  const double y_dbl = value_of(y_vec[n]);
87  if (y_dbl <= 0)
88  return LOG_ZERO;
89  }
90 
91  size_t N = max_size(y, alpha, beta);
92  agrad::OperandsAndPartials<T_y, T_shape, T_scale> operands_and_partials(y, alpha, beta);
93 
94  using boost::math::lgamma;
95  using stan::math::multiply_log;
96  using boost::math::digamma;
97 
98  DoubleVectorView<include_summand<propto,T_y,T_shape>::value,T_y>
99  log_y(length(y));
100  DoubleVectorView<include_summand<propto,T_y,T_scale>::value,T_y>
101  inv_y(length(y));
102  for(size_t n = 0; n < length(y); n++) {
103  if (include_summand<propto,T_y,T_shape>::value)
104  if (value_of(y_vec[n]) > 0)
105  log_y[n] = log(value_of(y_vec[n]));
106  if (include_summand<propto,T_y,T_scale>::value)
107  inv_y[n] = 1.0 / value_of(y_vec[n]);
108  }
109 
110  DoubleVectorView<include_summand<propto,T_shape>::value,T_shape>
111  lgamma_alpha(length(alpha));
112  DoubleVectorView<!is_constant_struct<T_shape>::value,T_shape>
113  digamma_alpha(length(alpha));
114  for (size_t n = 0; n < length(alpha); n++) {
115  if (include_summand<propto,T_shape>::value)
116  lgamma_alpha[n] = lgamma(value_of(alpha_vec[n]));
117  if (!is_constant_struct<T_shape>::value)
118  digamma_alpha[n] = digamma(value_of(alpha_vec[n]));
119  }
120 
121  DoubleVectorView<include_summand<propto,T_shape,T_scale>::value,T_scale>
122  log_beta(length(beta));
123  if (include_summand<propto,T_shape,T_scale>::value)
124  for (size_t n = 0; n < length(beta); n++)
125  log_beta[n] = log(value_of(beta_vec[n]));
126 
127  for (size_t n = 0; n < N; n++) {
128  // pull out values of arguments
129  const double alpha_dbl = value_of(alpha_vec[n]);
130  const double beta_dbl = value_of(beta_vec[n]);
131 
132  if (include_summand<propto,T_shape>::value)
133  logp -= lgamma_alpha[n];
134  if (include_summand<propto,T_shape,T_scale>::value)
135  logp += alpha_dbl * log_beta[n];
136  if (include_summand<propto,T_y,T_shape>::value)
137  logp -= (alpha_dbl+1.0) * log_y[n];
138  if (include_summand<propto,T_y,T_scale>::value)
139  logp -= beta_dbl * inv_y[n];
140 
141  // gradients
142  if (!is_constant<typename is_vector<T_y>::type>::value)
143  operands_and_partials.d_x1[n] += -(alpha_dbl+1) * inv_y[n] + beta_dbl * inv_y[n] * inv_y[n];
144  if (!is_constant<typename is_vector<T_shape>::type>::value)
145  operands_and_partials.d_x2[n] += -digamma_alpha[n] + log_beta[n] - log_y[n];
146  if (!is_constant<typename is_vector<T_scale>::type>::value)
147  operands_and_partials.d_x3[n] += alpha_dbl / beta_dbl - inv_y[n];
148  }
149  return operands_and_partials.to_var(logp);
150  }
151 
152  template <bool propto,
153  typename T_y, typename T_shape, typename T_scale>
154  inline
155  typename return_type<T_y,T_shape,T_scale>::type
156  inv_gamma_log(const T_y& y, const T_shape& alpha, const T_scale& beta) {
157  return inv_gamma_log<propto>(y,alpha,beta,stan::math::default_policy());
158  }
159 
160  template <typename T_y, typename T_shape, typename T_scale,
161  class Policy>
162  inline
163  typename return_type<T_y,T_shape,T_scale>::type
164  inv_gamma_log(const T_y& y, const T_shape& alpha, const T_scale& beta,
165  const Policy&) {
166  return inv_gamma_log<false>(y,alpha,beta,Policy());
167  }
168 
169  template <typename T_y, typename T_shape, typename T_scale>
170  inline
171  typename return_type<T_y,T_shape,T_scale>::type
172  inv_gamma_log(const T_y& y, const T_shape& alpha, const T_scale& beta) {
173  return inv_gamma_log<false>(y,alpha,beta,stan::math::default_policy());
174  }
175 
176 
177 
178  }
179 }
180 
181 #endif
     [ Stan Home Page ] © 2011–2012, Stan Development Team.