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>
6 #include <stan/math/prim/scal/fun/log_sum_exp.hpp>
40 void stable_sum(
double a1,
double log_w1,
double a2,
double log_w2,
41 double& sum_a,
double& log_sum_w) {
42 if (log_w2 > log_w1) {
43 double e = std::exp(log_w1 - log_w2);
44 sum_a = (e * a1 + a2) / (1 + e);
45 log_sum_w = log_w2 + std::log(1 + e);
47 double e = std::exp(log_w2 - log_w1);
48 sum_a = (a1 + e * a2) / (1 + e);
49 log_sum_w = log_w1 + std::log(1 + e);
57 template <
class Model,
template<
class,
class>
class Hamiltonian,
58 template<
class>
class Integrator,
class BaseRNG>
62 :
base_hmc<Model, Hamiltonian, Integrator, BaseRNG>(model, rng),
106 info_writer, error_writer);
107 double log_sum_weight = 0;
111 double sum_metro_prob = 1;
119 bool valid_subtree =
false;
120 double ave_subtree = 0;
121 double log_sum_weight_subtree
122 = -std::numeric_limits<double>::infinity();
125 this->
z_.ps_point::operator=(z_plus);
128 ave_subtree, log_sum_weight_subtree,
129 H0, 1, n_leapfrog, sum_metro_prob,
130 info_writer, error_writer);
131 z_plus.ps_point::operator=(this->
z_);
133 this->
z_.ps_point::operator=(z_minus);
136 ave_subtree, log_sum_weight_subtree,
137 H0, -1, n_leapfrog, sum_metro_prob,
138 info_writer, error_writer);
139 z_minus.ps_point::operator=(this->
z_);
142 if (!valid_subtree)
break;
144 ave_subtree, log_sum_weight_subtree,
145 ave, log_sum_weight);
151 = std::exp(log_sum_weight_subtree - log_sum_weight);
153 z_sample = z_propose;
165 = sum_metro_prob /
static_cast<double>(n_leapfrog + 1);
167 this->
z_.ps_point::operator=(z_sample);
169 return sample(this->
z_.q, -this->z_.V, accept_prob);
173 names.push_back(
"stepsize__");
174 names.push_back(
"treedepth__");
175 names.push_back(
"n_leapfrog__");
176 names.push_back(
"divergent__");
177 names.push_back(
"energy__");
182 values.push_back(this->
depth_);
185 values.push_back(this->
energy_);
205 double& ave,
double& log_sum_weight,
206 double H0,
double sign,
int& n_leapfrog,
207 double& sum_metro_prob,
214 info_writer, error_writer);
218 if (boost::math::isnan(h))
219 h = std::numeric_limits<double>::infinity();
229 ave, log_sum_weight);
234 sum_metro_prob += std::exp(H0 - h);
236 z_propose = this->
z_;
244 double log_sum_weight_left = -std::numeric_limits<double>::infinity();
248 ave_left, log_sum_weight_left,
249 H0, sign, n_leapfrog, sum_metro_prob,
250 info_writer, error_writer);
252 if (!valid_left)
return false;
254 ave_left, log_sum_weight_left,
255 ave, log_sum_weight);
259 double ave_right = 0;
260 double log_sum_weight_right = -std::numeric_limits<double>::infinity();
264 ave_right, log_sum_weight_right,
265 H0, sign, n_leapfrog, sum_metro_prob,
266 info_writer, error_writer);
268 if (!valid_right)
return false;
270 ave_right, log_sum_weight_right,
271 ave, log_sum_weight);
275 double log_sum_weight_subtree;
277 ave_right, log_sum_weight_right,
278 ave_subtree, log_sum_weight_subtree);
281 = std::exp(log_sum_weight_right - log_sum_weight_subtree);
283 z_propose = z_propose_right;
285 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.
double cont_params(int k) const
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)
Hamiltonian< Model, BaseRNG > hamiltonian_
Integrator< Hamiltonian< Model, BaseRNG > > integrator_