1 #ifndef STAN_MCMC_HMC_HAMILTONIANS_SOFTABS_METRIC_HPP 2 #define STAN_MCMC_HMC_HAMILTONIANS_SOFTABS_METRIC_HPP 4 #include <stan/math/mix/mat.hpp> 7 #include <boost/random/variate_generator.hpp> 8 #include <boost/random/normal_distribution.hpp> 13 template <
typename Model>
18 softabs_fun(
const Model& m, std::ostream* out): model_(m), o_(out) {}
21 T
operator()(Eigen::Matrix<T, Eigen::Dynamic, 1>& x)
const {
22 return model_.template log_prob<true, true, T>(x,
o_);
27 template <
class Model,
class BaseRNG>
31 typedef typename stan::math::index_type<Eigen::VectorXd>::type idx_t;
41 Eigen::VectorXd Qp = z.
eigen_deco.eigenvectors().transpose() * z.
p;
53 - z.
q.dot(dtau_dq(z, info_writer, error_writer)
54 + dphi_dq(z, info_writer, error_writer));
63 Eigen::MatrixXd A = a.asDiagonal()
65 Eigen::MatrixXd B = z.
pseudo_j.selfadjointView<Eigen::Lower>() * A;
66 Eigen::MatrixXd C = A.transpose() * B;
68 Eigen::VectorXd b(z.
q.size());
87 Eigen::MatrixXd A = a.asDiagonal()
89 Eigen::MatrixXd B = z.
eigen_deco.eigenvectors() * A;
94 return - 0.5 * a + z.
g;
98 boost::variate_generator<BaseRNG&, boost::normal_distribution<> >
99 rand_unit_gaus(rng, boost::normal_distribution<>());
101 Eigen::VectorXd a(z.
p.size());
103 for (idx_t n = 0; n < z.
p.size(); ++n)
113 update_metric(z, info_writer, error_writer);
114 update_metric_gradient(z, info_writer, error_writer);
122 stan::math::hessian<double>(
133 for (idx_t i = 0; i < z.
q.size(); ++i) {
134 double lambda = z.
eigen_deco.eigenvalues()(i);
135 double alpha_lambda = z.
alpha * lambda;
137 double softabs_lambda = 0;
141 if (std::fabs(alpha_lambda) < lower_softabs_thresh) {
142 softabs_lambda = (1.0
143 + (1.0 / 3.0) * alpha_lambda * alpha_lambda)
145 }
else if (std::fabs(alpha_lambda) > upper_softabs_thresh) {
146 softabs_lambda = std::fabs(lambda);
148 softabs_lambda = lambda / std::tanh(alpha_lambda);
157 for (idx_t i = 0; i < z.
q.size(); ++i)
166 for (idx_t i = 0; i < z.
q.size(); ++i) {
167 for (idx_t j = 0; j <= i; ++j) {
171 if (std::fabs(delta) < jacobian_thresh) {
172 double lambda = z.
eigen_deco.eigenvalues()(i);
173 double alpha_lambda = z.
alpha * lambda;
177 if (std::fabs(alpha_lambda) < lower_softabs_thresh) {
178 z.
pseudo_j(i, j) = (2.0 / 3.0) * alpha_lambda
179 * (1.0 - (2.0 / 15.0)
180 * alpha_lambda * alpha_lambda);
181 }
else if (std::fabs(alpha_lambda) > upper_softabs_thresh) {
182 z.
pseudo_j(i, j) = lambda > 0 ? 1 : -1;
184 double sdx = std::sinh(alpha_lambda) / lambda;
186 - z.
alpha / (sdx * sdx) ) / lambda;
200 update_metric_gradient(z, info_writer, error_writer);
217 template <
class Model,
class BaseRNG>
220 template <
class Model,
class BaseRNG>
223 template <
class Model,
class BaseRNG>
Eigen::VectorXd dphi_dq(softabs_point &z, interface_callbacks::writer::base_writer &info_writer, interface_callbacks::writer::base_writer &error_writer)
Probability, optimization and sampling library.
Eigen::VectorXd dtau_dq(softabs_point &z, interface_callbacks::writer::base_writer &info_writer, interface_callbacks::writer::base_writer &error_writer)
Eigen::VectorXd softabs_lambda_inv
static double jacobian_thresh
Eigen::VectorXd softabs_lambda
static double lower_softabs_thresh
void init(softabs_point &z, interface_callbacks::writer::base_writer &info_writer, interface_callbacks::writer::base_writer &error_writer)
double T(softabs_point &z)
double dG_dt(softabs_point &z, interface_callbacks::writer::base_writer &info_writer, interface_callbacks::writer::base_writer &error_writer)
void update_metric(softabs_point &z, interface_callbacks::writer::base_writer &info_writer, interface_callbacks::writer::base_writer &error_writer)
Eigen::SelfAdjointEigenSolver< Eigen::MatrixXd > eigen_deco
double tau(softabs_point &z)
Point in a phase space with a base Riemannian manifold with SoftAbs metric.
base_writer is an abstract base class defining the interface for Stan writer callbacks.
T operator()(Eigen::Matrix< T, Eigen::Dynamic, 1 > &x) const
void update_gradients(softabs_point &z, interface_callbacks::writer::base_writer &info_writer, interface_callbacks::writer::base_writer &error_writer)
void update_metric_gradient(softabs_point &z, interface_callbacks::writer::base_writer &info_writer, interface_callbacks::writer::base_writer &error_writer)
void grad_tr_mat_times_hessian(const M &model, const Eigen::Matrix< double, Eigen::Dynamic, 1 > &x, const Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > &X, Eigen::Matrix< double, Eigen::Dynamic, 1 > &grad_tr_X_hess_f, std::ostream *msgs=0)
Eigen::VectorXd dtau_dp(softabs_point &z)
softabs_fun(const Model &m, std::ostream *out)
double phi(softabs_point &z)
static double upper_softabs_thresh
softabs_metric(const Model &model)
void sample_p(softabs_point &z, BaseRNG &rng)