1 #ifndef STAN_MCMC_HMC_NUTS_BASE_XHMC_HPP 2 #define STAN_MCMC_HMC_NUTS_BASE_XHMC_HPP 5 #include <boost/math/special_functions/fpclassify.hpp> 39 void stable_sum(
double a1,
double log_w1,
double a2,
double log_w2,
40 double& sum_a,
double& log_sum_w) {
41 if (log_w2 > log_w1) {
42 double e = std::exp(log_w1 - log_w2);
43 sum_a = (e * a1 + a2) / (1 + e);
44 log_sum_w = log_w2 + std::log(1 + e);
46 double e = std::exp(log_w2 - log_w1);
47 sum_a = (a1 + e * a2) / (1 + e);
48 log_sum_w = log_w1 + std::log(1 + e);
56 template <
class Model,
template<
class,
class>
class Hamiltonian,
57 template<
class>
class Integrator,
class BaseRNG>
61 :
base_hmc<Model, Hamiltonian, Integrator, BaseRNG>(model, rng),
105 info_writer, error_writer);
106 double log_sum_weight = 0;
110 double sum_metro_prob = 1;
118 bool valid_subtree =
false;
119 double ave_subtree = 0;
120 double log_sum_weight_subtree
121 = -std::numeric_limits<double>::infinity();
124 this->
z_.ps_point::operator=(z_plus);
127 ave_subtree, log_sum_weight_subtree,
128 H0, 1, n_leapfrog, sum_metro_prob,
129 info_writer, error_writer);
130 z_plus.ps_point::operator=(this->
z_);
132 this->
z_.ps_point::operator=(z_minus);
135 ave_subtree, log_sum_weight_subtree,
136 H0, -1, n_leapfrog, sum_metro_prob,
137 info_writer, error_writer);
138 z_minus.ps_point::operator=(this->
z_);
141 if (!valid_subtree)
break;
143 ave_subtree, log_sum_weight_subtree,
144 ave, log_sum_weight);
150 = std::exp(log_sum_weight_subtree - log_sum_weight);
152 z_sample = z_propose;
164 = sum_metro_prob /
static_cast<double>(n_leapfrog + 1);
166 this->
z_.ps_point::operator=(z_sample);
168 return sample(this->
z_.q, -this->z_.V, accept_prob);
172 names.push_back(
"stepsize__");
173 names.push_back(
"treedepth__");
174 names.push_back(
"n_leapfrog__");
175 names.push_back(
"divergent__");
176 names.push_back(
"energy__");
181 values.push_back(this->
depth_);
184 values.push_back(this->
energy_);
204 double& ave,
double& log_sum_weight,
205 double H0,
double sign,
int& n_leapfrog,
206 double& sum_metro_prob,
213 info_writer, error_writer);
217 if (boost::math::isnan(h))
218 h = std::numeric_limits<double>::infinity();
228 ave, log_sum_weight);
233 sum_metro_prob += std::exp(H0 - h);
235 z_propose = this->
z_;
243 double log_sum_weight_left = -std::numeric_limits<double>::infinity();
247 ave_left, log_sum_weight_left,
248 H0, sign, n_leapfrog, sum_metro_prob,
249 info_writer, error_writer);
251 if (!valid_left)
return false;
253 ave_left, log_sum_weight_left,
254 ave, log_sum_weight);
258 double ave_right = 0;
259 double log_sum_weight_right = -std::numeric_limits<double>::infinity();
263 ave_right, log_sum_weight_right,
264 H0, sign, n_leapfrog, sum_metro_prob,
265 info_writer, error_writer);
267 if (!valid_right)
return false;
269 ave_right, log_sum_weight_right,
270 ave, log_sum_weight);
274 double log_sum_weight_subtree;
276 ave_right, log_sum_weight_right,
277 ave_subtree, log_sum_weight_subtree);
280 = std::exp(log_sum_weight_right - log_sum_weight_subtree);
282 z_propose = z_propose_right;
284 return std::fabs(ave_subtree) >=
x_delta_;
Hamiltonian< Model, BaseRNG >::PointType z_
base_xhmc(const Model &model, BaseRNG &rng)
sample transition(sample &init_sample, interface_callbacks::writer::base_writer &info_writer, interface_callbacks::writer::base_writer &error_writer)
Exhaustive Hamiltonian Monte Carlo (XHMC) with multinomial sampling.
void sample(stan::mcmc::base_mcmc *sampler, int num_warmup, int num_samples, int num_thin, int refresh, bool save, stan::services::sample::mcmc_writer< Model, SampleRecorder, DiagnosticRecorder, MessageRecorder > &mcmc_writer, stan::mcmc::sample &init_s, Model &model, RNG &base_rng, const std::string &prefix, const std::string &suffix, std::ostream &o, StartTransitionCallback &callback, interface_callbacks::writer::base_writer &info_writer, interface_callbacks::writer::base_writer &error_writer)
Probability, optimization and sampling library.
Point in a generic phase space.
void set_max_depth(int d)
int build_tree(int depth, ps_point &z_propose, double &ave, double &log_sum_weight, double H0, double sign, int &n_leapfrog, double &sum_metro_prob, interface_callbacks::writer::base_writer &info_writer, interface_callbacks::writer::base_writer &error_writer)
Recursively build a new subtree to completion or until the subtree becomes invalid.
base_writer is an abstract base class defining the interface for Stan writer callbacks.
void stable_sum(double a1, double log_w1, double a2, double log_w2, double &sum_a, double &log_sum_w)
a1 and a2 are running averages of the form and the weights are the respective normalizing constants...
void set_x_delta(double d)
void get_sampler_params(std::vector< double > &values)
void seed(const Eigen::VectorXd &q)
void get_sampler_param_names(std::vector< std::string > &names)
boost::uniform_01< BaseRNG & > rand_uniform_
void set_max_deltaH(double d)
double cont_params(int k) const
Hamiltonian< Model, BaseRNG > hamiltonian_
Integrator< Hamiltonian< Model, BaseRNG > > integrator_