1 #ifndef STAN_MCMC_HMC_UNIFORM_BASE_STATIC_UNIFORM_HPP 2 #define STAN_MCMC_HMC_UNIFORM_BASE_STATIC_UNIFORM_HPP 3 4 #include <stan/callbacks/logger.hpp> 5 #include <stan/mcmc/hmc/base_hmc.hpp> 6 #include <stan/mcmc/hmc/hamiltonians/ps_point.hpp> 7 #include <boost/random/uniform_int_distribution.hpp> 8 #include <cmath> 9 #include <limits> 10 #include <string> 11 #include <vector> 12 13 namespace stan { 14 15 namespace mcmc { 16 /** 17 * Hamiltonian Monte Carlo implementation that uniformly samples 18 * from trajectories with a static integration time 19 */ 20 template <class Model, template <class, class> class Hamiltonian, 21 template <class> class Integrator, class BaseRNG> 22 class base_static_uniform 23 : public base_hmc<Model, Hamiltonian, Integrator, BaseRNG> { 24 public: base_static_uniform(const Model & model,BaseRNG & rng)25 base_static_uniform(const Model& model, BaseRNG& rng) 26 : base_hmc<Model, Hamiltonian, Integrator, BaseRNG>(model, rng), 27 T_(1), 28 energy_(0) { 29 update_L_(); 30 } 31 ~base_static_uniform()32 ~base_static_uniform() {} 33 transition(sample & init_sample,callbacks::logger & logger)34 sample transition(sample& init_sample, callbacks::logger& logger) { 35 this->sample_stepsize(); 36 37 this->seed(init_sample.cont_params()); 38 39 this->hamiltonian_.sample_p(this->z_, this->rand_int_); 40 this->hamiltonian_.init(this->z_, logger); 41 42 ps_point z_init(this->z_); 43 double H0 = this->hamiltonian_.H(this->z_); 44 45 ps_point z_sample(this->z_); 46 double sum_prob = 1; 47 double sum_metro_prob = 1; 48 49 boost::random::uniform_int_distribution<> uniform(0, L_ - 1); 50 int Lp = uniform(this->rand_int_); 51 52 for (int l = 0; l < Lp; ++l) { 53 this->integrator_.evolve(this->z_, this->hamiltonian_, -this->epsilon_, 54 logger); 55 56 double h = this->hamiltonian_.H(this->z_); 57 if (std::isnan(h)) 58 h = std::numeric_limits<double>::infinity(); 59 60 double prob = std::exp(H0 - h); 61 sum_prob += prob; 62 sum_metro_prob += prob > 1 ? 1 : prob; 63 64 if (this->rand_uniform_() < prob / sum_prob) 65 z_sample = this->z_; 66 } 67 68 this->z_.ps_point::operator=(z_init); 69 70 for (int l = 0; l < L_ - 1 - Lp; ++l) { 71 this->integrator_.evolve(this->z_, this->hamiltonian_, this->epsilon_, 72 logger); 73 74 double h = this->hamiltonian_.H(this->z_); 75 if (std::isnan(h)) 76 h = std::numeric_limits<double>::infinity(); 77 78 double prob = std::exp(H0 - h); 79 sum_prob += prob; 80 sum_metro_prob += prob > 1 ? 1 : prob; 81 82 if (this->rand_uniform_() < prob / sum_prob) 83 z_sample = this->z_; 84 } 85 86 double accept_prob = sum_metro_prob / static_cast<double>(L_); 87 88 this->z_.ps_point::operator=(z_sample); 89 this->energy_ = this->hamiltonian_.H(this->z_); 90 return sample(this->z_.q, -this->hamiltonian_.V(this->z_), accept_prob); 91 } 92 get_sampler_param_names(std::vector<std::string> & names)93 void get_sampler_param_names(std::vector<std::string>& names) { 94 names.push_back("stepsize__"); 95 names.push_back("int_time__"); 96 names.push_back("energy__"); 97 } 98 get_sampler_params(std::vector<double> & values)99 void get_sampler_params(std::vector<double>& values) { 100 values.push_back(this->epsilon_); 101 values.push_back(this->T_); 102 values.push_back(this->energy_); 103 } 104 set_nominal_stepsize_and_T(const double e,const double t)105 void set_nominal_stepsize_and_T(const double e, const double t) { 106 if (e > 0 && t > 0) { 107 this->nom_epsilon_ = e; 108 T_ = t; 109 update_L_(); 110 } 111 } 112 set_nominal_stepsize_and_L(const double e,const int l)113 void set_nominal_stepsize_and_L(const double e, const int l) { 114 if (e > 0 && l > 0) { 115 this->nom_epsilon_ = e; 116 L_ = l; 117 T_ = this->nom_epsilon_ * L_; 118 } 119 } 120 set_T(const double t)121 void set_T(const double t) { 122 if (t > 0) { 123 T_ = t; 124 update_L_(); 125 } 126 } 127 set_nominal_stepsize(const double e)128 void set_nominal_stepsize(const double e) { 129 if (e > 0) { 130 this->nom_epsilon_ = e; 131 update_L_(); 132 } 133 } 134 get_T()135 double get_T() { return this->T_; } 136 get_L()137 int get_L() { return this->L_; } 138 139 protected: 140 double T_; 141 int L_; 142 double energy_; 143 update_L_()144 void update_L_() { 145 L_ = static_cast<int>(T_ / this->nom_epsilon_); 146 L_ = L_ < 1 ? 1 : L_; 147 } 148 }; 149 } // namespace mcmc 150 } // namespace stan 151 #endif 152