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>
9 #include <boost/random/variate_generator.hpp>
10 #include <boost/random/normal_distribution.hpp>
15 template <
typename Model>
20 softabs_fun(
const Model& m, std::ostream* out): model_(m), o_(out) {}
23 T
operator()(Eigen::Matrix<T, Eigen::Dynamic, 1>& x)
const {
24 return model_.template log_prob<true, true, T>(x,
o_);
29 template <
class Model,
class BaseRNG>
33 typedef typename stan::math::index_type<Eigen::VectorXd>::type idx_t;
43 Eigen::VectorXd Qp = z.
eigen_deco.eigenvectors().transpose() * z.
p;
55 - z.
q.dot(
dtau_dq(z, info_writer, error_writer)
56 +
dphi_dq(z, info_writer, error_writer));
65 Eigen::MatrixXd A = a.asDiagonal()
67 Eigen::MatrixXd B = z.
pseudo_j.selfadjointView<Eigen::Lower>() * A;
68 Eigen::MatrixXd C = A.transpose() * B;
70 Eigen::VectorXd b(z.
q.size());
89 Eigen::MatrixXd A = a.asDiagonal()
91 Eigen::MatrixXd B = z.
eigen_deco.eigenvectors() * A;
96 return - 0.5 * a + z.
g;
100 boost::variate_generator<BaseRNG&, boost::normal_distribution<> >
101 rand_unit_gaus(rng, boost::normal_distribution<>());
103 Eigen::VectorXd a(z.
p.size());
105 for (idx_t n = 0; n < z.
p.size(); ++n)
124 stan::math::hessian<double>(
135 for (idx_t i = 0; i < z.
q.size(); ++i) {
136 double lambda = z.
eigen_deco.eigenvalues()(i);
137 double alpha_lambda = z.
alpha * lambda;
139 double softabs_lambda = 0;
144 softabs_lambda = (1.0
145 + (1.0 / 3.0) * alpha_lambda * alpha_lambda)
148 softabs_lambda = std::fabs(lambda);
150 softabs_lambda = lambda / std::tanh(alpha_lambda);
159 for (idx_t i = 0; i < z.
q.size(); ++i)
168 for (idx_t i = 0; i < z.
q.size(); ++i) {
169 for (idx_t j = 0; j <= i; ++j) {
174 double lambda = z.
eigen_deco.eigenvalues()(i);
175 double alpha_lambda = z.
alpha * lambda;
180 z.
pseudo_j(i, j) = (2.0 / 3.0) * alpha_lambda
181 * (1.0 - (2.0 / 15.0)
182 * alpha_lambda * alpha_lambda);
184 z.
pseudo_j(i, j) = lambda > 0 ? 1 : -1;
186 double sdx = std::sinh(alpha_lambda) / lambda;
188 - z.
alpha / (sdx * sdx) ) / lambda;
219 template <
class Model,
class BaseRNG>
222 template <
class Model,
class BaseRNG>
225 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.
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)
double V(softabs_point &z)
T operator()(Eigen::Matrix< T, Eigen::Dynamic, 1 > &x) const
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)