Automatic Differentiation
 
Loading...
Searching...
No Matches
log_mix.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_FUN_LOG_MIX_HPP
2#define STAN_MATH_PRIM_FUN_LOG_MIX_HPP
3
16#include <vector>
17#include <cmath>
18
19namespace stan {
20namespace math {
21
39template <typename T_theta, typename T_lambda1, typename T_lambda2,
40 require_all_arithmetic_t<T_theta, T_lambda1, T_lambda2>* = nullptr>
41inline double log_mix(T_theta theta, T_lambda1 lambda1, T_lambda2 lambda2) {
42 using std::log;
43 check_not_nan("log_mix", "lambda1", lambda1);
44 check_not_nan("log_mix", "lambda2", lambda2);
45 check_bounded("log_mix", "theta", theta, 0, 1);
46 return log_sum_exp(log(theta) + lambda1, log1m(theta) + lambda2);
47}
48
75template <typename T_theta, typename T_lam,
78 const T_lam& lambda) {
79 static constexpr const char* function = "log_mix";
80 using T_partials_return = partials_return_t<T_theta, T_lam>;
81 using T_partials_vec =
82 typename Eigen::Matrix<T_partials_return, Eigen::Dynamic, 1>;
83 using T_theta_ref = ref_type_t<T_theta>;
84 using T_lam_ref = ref_type_t<T_lam>;
85
86 check_consistent_sizes(function, "theta", theta, "lambda", lambda);
87 T_theta_ref theta_ref = theta;
88 T_lam_ref lambda_ref = lambda;
89 check_bounded(function, "theta", theta_ref, 0, 1);
90 check_finite(function, "lambda", lambda_ref);
91
92 const auto& theta_dbl
94 const auto& lam_dbl
96
97 T_partials_return logp = log_sum_exp(log(theta_dbl) + lam_dbl);
98
99 auto ops_partials = make_partials_propagator(theta_ref, lambda_ref);
101 T_partials_vec theta_deriv = (lam_dbl.array() - logp).exp();
103 partials<1>(ops_partials) = theta_deriv.cwiseProduct(theta_dbl);
104 }
106 partials<0>(ops_partials) = std::move(theta_deriv);
107 }
108 }
109 return ops_partials.build(logp);
110}
111
143template <typename T_theta, typename T_lam, require_vector_t<T_lam>* = nullptr>
145 const T_theta& theta, const std::vector<T_lam>& lambda) {
146 static constexpr const char* function = "log_mix";
147 using T_partials_return = partials_return_t<T_theta, std::vector<T_lam>>;
148 using T_partials_vec =
149 typename Eigen::Matrix<T_partials_return, Eigen::Dynamic, 1>;
150 using T_partials_mat =
151 typename Eigen::Matrix<T_partials_return, Eigen::Dynamic, Eigen::Dynamic>;
152 using T_theta_ref = ref_type_t<T_theta>;
153
154 const int N = stan::math::size(lambda);
155 const int M = theta.size();
156
157 T_theta_ref theta_ref = theta;
158 check_bounded(function, "theta", theta_ref, 0, 1);
159 for (int n = 0; n < N; ++n) {
160 check_not_nan(function, "lambda", lambda[n]);
161 check_finite(function, "lambda", lambda[n]);
162 check_consistent_sizes(function, "theta", theta, "lambda", lambda[n]);
163 }
164
165 const auto& theta_dbl
167
168 T_partials_mat lam_dbl(M, N);
169 for (int n = 0; n < N; ++n) {
170 lam_dbl.col(n) = value_of(as_column_vector_or_scalar(lambda[n]));
171 }
172
173 T_partials_mat logp_tmp = lam_dbl.colwise() + log(theta_dbl);
174 T_partials_vec logp(N);
175 for (int n = 0; n < N; ++n) {
176 logp[n] = log_sum_exp(logp_tmp.col(n));
177 }
178
179 auto ops_partials = make_partials_propagator(theta_ref, lambda);
181 T_partials_mat derivs = exp(lam_dbl.rowwise() - logp.transpose());
183 partials<0>(ops_partials) = derivs.rowwise().sum();
184 }
186 for (int n = 0; n < N; ++n) {
187 as_column_vector_or_scalar(partials_vec<1>(ops_partials)[n])
188 = derivs.col(n).cwiseProduct(theta_dbl);
189 }
190 }
191 }
192 return ops_partials.build(logp.sum());
193}
194
195} // namespace math
196} // namespace stan
197#endif
auto as_column_vector_or_scalar(T &&a)
as_column_vector_or_scalar of a kernel generator expression.
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
require_any_t< is_vector< std::decay_t< Types > >... > require_any_vector_t
Require any of the types satisfy is_vector.
void check_bounded(const char *function, const char *name, const T_y &y, const T_low &low, const T_high &high)
Check if the value is between the low and high values, inclusively.
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition value_of.hpp:18
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 > log_mix(const fvar< T > &theta, const fvar< T > &lambda1, const fvar< T > &lambda2)
Return the log mixture density with specified mixing proportion and log densities and its derivative ...
Definition log_mix.hpp:98
ref_type_t< T && > to_ref(T &&a)
This evaluates expensive Eigen expressions.
Definition to_ref.hpp:17
void check_finite(const char *function, const char *name, const T_y &y)
Return true if all values in y are finite.
void check_not_nan(const char *function, const char *name, const T_y &y)
Check if y is not NaN.
fvar< T > log1m(const fvar< T > &x)
Definition log1m.hpp:12
auto make_partials_propagator(Ops &&... ops)
Construct an partials_propagator.
fvar< T > log_sum_exp(const fvar< T > &x1, const fvar< T > &x2)
fvar< T > exp(const fvar< T > &x)
Definition exp.hpp:13
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...