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