Automatic Differentiation
 
Loading...
Searching...
No Matches
scaled_inv_chi_square_lpdf.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_PROB_SCALED_INV_CHI_SQUARE_LPDF_HPP
2#define STAN_MATH_PRIM_PROB_SCALED_INV_CHI_SQUARE_LPDF_HPP
3
17#include <cmath>
18
19namespace stan {
20namespace math {
21
44template <bool propto, typename T_y, typename T_dof, typename T_scale,
46 T_y, T_dof, T_scale>* = nullptr>
48 const T_y& y, const T_dof& nu, const T_scale& s) {
49 using T_partials_return = partials_return_t<T_y, T_dof, T_scale>;
50 using std::log;
51 using T_y_ref = ref_type_t<T_y>;
52 using T_nu_ref = ref_type_t<T_dof>;
53 using T_s_ref = ref_type_t<T_scale>;
54 static constexpr const char* function = "scaled_inv_chi_square_lpdf";
55 check_consistent_sizes(function, "Random variable", y,
56 "Degrees of freedom parameter", nu, "Scale parameter",
57 s);
58 T_y_ref y_ref = y;
59 T_nu_ref nu_ref = nu;
60 T_s_ref s_ref = s;
61 check_not_nan(function, "Random variable", y_ref);
62 check_positive_finite(function, "Degrees of freedom parameter", nu_ref);
63 check_positive_finite(function, "Scale parameter", s_ref);
64
65 if (size_zero(y, nu, s)) {
66 return 0;
67 }
69 return 0;
70 }
71
72 T_partials_return logp(0);
73 auto ops_partials = make_partials_propagator(y_ref, nu_ref, s_ref);
74
75 scalar_seq_view<T_y_ref> y_vec(y_ref);
76 scalar_seq_view<T_nu_ref> nu_vec(nu_ref);
77 scalar_seq_view<T_s_ref> s_vec(s_ref);
78 size_t N = max_size(y, nu, s);
79
80 for (size_t n = 0; n < N; n++) {
81 if (y_vec.val(n) <= 0) {
82 return LOG_ZERO;
83 }
84 }
85
87 T_partials_return, T_dof>
88 half_nu(math::size(nu));
89 for (size_t i = 0; i < stan::math::size(nu); i++) {
91 half_nu[i] = 0.5 * nu_vec.val(i);
92 }
93 }
94
96 T_y>
97 log_y(math::size(y));
98 for (size_t i = 0; i < stan::math::size(y); i++) {
100 log_y[i] = log(y_vec.val(i));
101 }
102 }
103
105 T_partials_return, T_y>
106 inv_y(math::size(y));
107 for (size_t i = 0; i < stan::math::size(y); i++) {
109 inv_y[i] = 1.0 / y_vec.val(i);
110 }
111 }
112
114 T_partials_return, T_scale>
115 log_s(math::size(s));
116 for (size_t i = 0; i < stan::math::size(s); i++) {
118 log_s[i] = log(s_vec.val(i));
119 }
120 }
121
123 log_half_nu(math::size(nu));
125 lgamma_half_nu(math::size(nu));
126 VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
127 digamma_half_nu_over_two(math::size(nu));
128 for (size_t i = 0; i < stan::math::size(nu); i++) {
130 lgamma_half_nu[i] = lgamma(half_nu[i]);
131 }
133 log_half_nu[i] = log(half_nu[i]);
134 }
136 digamma_half_nu_over_two[i] = digamma(half_nu[i]) * 0.5;
137 }
138 }
139
140 for (size_t n = 0; n < N; n++) {
141 const T_partials_return s_dbl = s_vec.val(n);
142 const T_partials_return nu_dbl = nu_vec.val(n);
144 logp += half_nu[n] * log_half_nu[n] - lgamma_half_nu[n];
145 }
147 logp += nu_dbl * log_s[n];
148 }
150 logp -= (half_nu[n] + 1.0) * log_y[n];
151 }
153 logp -= half_nu[n] * s_dbl * s_dbl * inv_y[n];
154 }
155
157 partials<0>(ops_partials)[n]
158 += -(half_nu[n] + 1.0) * inv_y[n]
159 + half_nu[n] * s_dbl * s_dbl * inv_y[n] * inv_y[n];
160 }
162 partials<1>(ops_partials)[n]
163 += 0.5 * log_half_nu[n] + 0.5 - digamma_half_nu_over_two[n] + log_s[n]
164 - 0.5 * log_y[n] - 0.5 * s_dbl * s_dbl * inv_y[n];
165 }
167 partials<2>(ops_partials)[n]
168 += nu_dbl / s_dbl - nu_dbl * inv_y[n] * s_dbl;
169 }
170 }
171 return ops_partials.build(logp);
172}
173
174template <typename T_y, typename T_dof, typename T_scale>
176 const T_y& y, const T_dof& nu, const T_scale& s) {
177 return scaled_inv_chi_square_lpdf<false>(y, nu, s);
178}
179
180} // namespace math
181} // namespace stan
182#endif
VectorBuilder allocates type T1 values to be used as intermediate values.
scalar_seq_view provides a uniform sequence-like wrapper around either a scalar or a sequence of scal...
require_all_not_t< is_nonscalar_prim_or_rev_kernel_expression< std::decay_t< Types > >... > require_all_not_nonscalar_prim_or_rev_kernel_expression_t
Require none of the types satisfy is_nonscalar_prim_or_rev_kernel_expression.
return_type_t< T_y_cl, T_dof_cl, T_scale_cl > scaled_inv_chi_square_lpdf(const T_y_cl &y, const T_dof_cl &nu, const T_scale_cl &s)
The log of a scaled inverse chi-squared density for y with the specified degrees of freedom parameter...
typename return_type< Ts... >::type return_type_t
Convenience type for the return type of the specified template parameters.
int64_t size(const T &m)
Returns the size (number of the elements) of a matrix_cl or var_value<matrix_cl<T>>.
Definition size.hpp:19
static constexpr double LOG_ZERO
The natural logarithm of 0, .
Definition constants.hpp:68
bool size_zero(const T &x)
Returns 1 if input is of length 0, returns 0 otherwise.
Definition size_zero.hpp:19
fvar< T > log(const fvar< T > &x)
Definition log.hpp:15
void check_consistent_sizes(const char *)
Trivial no input case, this function is a no-op.
fvar< T > lgamma(const fvar< T > &x)
Return the natural logarithm of the gamma function applied to the specified argument.
Definition lgamma.hpp:21
void check_not_nan(const char *function, const char *name, const T_y &y)
Check if y is not NaN.
int64_t max_size(const T1 &x1, const Ts &... xs)
Calculate the size of the largest input.
Definition max_size.hpp:20
auto make_partials_propagator(Ops &&... ops)
Construct an partials_propagator.
void check_positive_finite(const char *function, const char *name, const T_y &y)
Check if y is positive and finite.
fvar< T > digamma(const fvar< T > &x)
Return the derivative of the log gamma function at the specified argument.
Definition digamma.hpp:23
typename ref_type_if< true, T >::type ref_type_t
Definition ref_type.hpp:55
typename partials_return_type< Args... >::type partials_return_t
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...
Extends std::true_type when instantiated with zero or more template parameters, all of which extend t...
Template metaprogram to calculate whether a summand needs to be included in a proportional (log) prob...