base_hmc.hpp
Go to the documentation of this file.
1 #ifndef STAN_MCMC_HMC_BASE_HMC_HPP
2 #define STAN_MCMC_HMC_BASE_HMC_HPP
3 
8 #include <boost/math/special_functions/fpclassify.hpp>
9 #include <boost/random/variate_generator.hpp>
10 #include <boost/random/uniform_01.hpp>
11 #include <cmath>
12 #include <limits>
13 #include <stdexcept>
14 #include <string>
15 #include <vector>
16 
17 namespace stan {
18  namespace mcmc {
19 
20  template <class Model,
21  template<class, class> class Hamiltonian,
22  template<class> class Integrator,
23  class BaseRNG>
24  class base_hmc : public base_mcmc {
25  public:
26  base_hmc(const Model &model, BaseRNG& rng)
27  : base_mcmc(),
28  z_(model.num_params_r()),
29  integrator_(),
30  hamiltonian_(model),
31  rand_int_(rng),
33  nom_epsilon_(0.1),
35  epsilon_jitter_(0.0) {}
36 
37  /**
38  * format and write stepsize
39  */
40  void
42  std::stringstream nominal_stepsize;
43  nominal_stepsize << "Step size = " << get_nominal_stepsize();
44  writer(nominal_stepsize.str());
45  }
46 
47  /**
48  * write elements of mass matrix
49  */
50  void
52  z_.write_metric(writer);
53  }
54 
55  /**
56  * write stepsize and elements of mass matrix
57  */
58  void
60  write_sampler_stepsize(writer);
61  write_sampler_metric(writer);
62  }
63 
64  void get_sampler_diagnostic_names(std::vector<std::string>& model_names,
65  std::vector<std::string>& names) {
66  z_.get_param_names(model_names, names);
67  }
68 
69  void get_sampler_diagnostics(std::vector<double>& values) {
70  z_.get_params(values);
71  }
72 
73  void seed(const Eigen::VectorXd& q) {
74  z_.q = q;
75  }
76 
77  void
79  this->hamiltonian_.init(this->z_, logger);
80  }
81 
82  void
84  ps_point z_init(this->z_);
85 
86  // Skip initialization for extreme step sizes
87  if (this->nom_epsilon_ == 0 || this->nom_epsilon_ > 1e7)
88  return;
89 
90  this->hamiltonian_.sample_p(this->z_, this->rand_int_);
91  this->hamiltonian_.init(this->z_, logger);
92 
93  // Guaranteed to be finite if randomly initialized
94  double H0 = this->hamiltonian_.H(this->z_);
95 
96  this->integrator_.evolve(this->z_, this->hamiltonian_,
97  this->nom_epsilon_,
98  logger);
99 
100  double h = this->hamiltonian_.H(this->z_);
101  if (boost::math::isnan(h))
102  h = std::numeric_limits<double>::infinity();
103 
104  double delta_H = H0 - h;
105 
106  int direction = delta_H > std::log(0.8) ? 1 : -1;
107 
108  while (1) {
109  this->z_.ps_point::operator=(z_init);
110 
111  this->hamiltonian_.sample_p(this->z_, this->rand_int_);
112  this->hamiltonian_.init(this->z_, logger);
113 
114  double H0 = this->hamiltonian_.H(this->z_);
115 
116  this->integrator_.evolve(this->z_, this->hamiltonian_,
117  this->nom_epsilon_,
118  logger);
119 
120  double h = this->hamiltonian_.H(this->z_);
121  if (boost::math::isnan(h))
122  h = std::numeric_limits<double>::infinity();
123 
124  double delta_H = H0 - h;
125 
126  if ((direction == 1) && !(delta_H > std::log(0.8)))
127  break;
128  else if ((direction == -1) && !(delta_H < std::log(0.8)))
129  break;
130  else
131  this->nom_epsilon_
132  = direction == 1
133  ? 2.0 * this->nom_epsilon_
134  : 0.5 * this->nom_epsilon_;
135 
136  if (this->nom_epsilon_ > 1e7)
137  throw std::runtime_error("Posterior is improper. "
138  "Please check your model.");
139  if (this->nom_epsilon_ == 0)
140  throw std::runtime_error("No acceptably small step size could "
141  "be found. Perhaps the posterior is "
142  "not continuous?");
143  }
144 
145  this->z_.ps_point::operator=(z_init);
146  }
147 
148  typename Hamiltonian<Model, BaseRNG>::PointType& z() {
149  return z_;
150  }
151 
152  virtual void set_nominal_stepsize(double e) {
153  if (e > 0)
154  nom_epsilon_ = e;
155  }
156 
158  return this->nom_epsilon_;
159  }
160 
162  return this->epsilon_;
163  }
164 
165  virtual void set_stepsize_jitter(double j) {
166  if (j > 0 && j < 1)
167  epsilon_jitter_ = j;
168  }
169 
171  return this->epsilon_jitter_;
172  }
173 
175  this->epsilon_ = this->nom_epsilon_;
176  if (this->epsilon_jitter_)
177  this->epsilon_ *= 1.0
178  + this->epsilon_jitter_ * (2.0 * this->rand_uniform_() - 1.0);
179  }
180 
181  protected:
182  typename Hamiltonian<Model, BaseRNG>::PointType z_;
183  Integrator<Hamiltonian<Model, BaseRNG> > integrator_;
184  Hamiltonian<Model, BaseRNG> hamiltonian_;
185 
186  BaseRNG& rand_int_;
187 
188  // Uniform(0, 1) RNG
189  boost::uniform_01<BaseRNG&> rand_uniform_;
190 
191  double nom_epsilon_;
192  double epsilon_;
194  };
195 
196  } // mcmc
197 } // stan
198 #endif
void write_sampler_state(callbacks::writer &writer)
Definition: base_hmc.hpp:59
Hamiltonian< Model, BaseRNG >::PointType z_
Definition: base_hmc.hpp:182
void get_sampler_diagnostics(std::vector< double > &values)
Definition: base_hmc.hpp:69
void init_stepsize(callbacks::logger &logger)
Definition: base_hmc.hpp:83
void write_sampler_metric(callbacks::writer &writer)
Definition: base_hmc.hpp:51
virtual void set_nominal_stepsize(double e)
Definition: base_hmc.hpp:152
void write_sampler_stepsize(callbacks::writer &writer)
Definition: base_hmc.hpp:41
void get_sampler_diagnostic_names(std::vector< std::string > &model_names, std::vector< std::string > &names)
Definition: base_hmc.hpp:64
rosenbrock_model_namespace::rosenbrock_model Model
virtual void set_stepsize_jitter(double j)
Definition: base_hmc.hpp:165
bool isnan(const stan::math::var &v)
Definition: boost_isnan.hpp:20
double get_current_stepsize()
Definition: base_hmc.hpp:161
const double j
Definition: BetheBloch.cxx:29
void seed(const Eigen::VectorXd &q)
Definition: base_hmc.hpp:73
base_hmc(const Model &model, BaseRNG &rng)
Definition: base_hmc.hpp:26
double get_stepsize_jitter()
Definition: base_hmc.hpp:170
void init_hamiltonian(callbacks::logger &logger)
Definition: base_hmc.hpp:78
double get_nominal_stepsize()
Definition: base_hmc.hpp:157
boost::uniform_01< BaseRNG & > rand_uniform_
Definition: base_hmc.hpp:189
Hamiltonian< Model, BaseRNG >::PointType & z()
Definition: base_hmc.hpp:148
Float_t e
Definition: plot.C:35
Hamiltonian< Model, BaseRNG > hamiltonian_
Definition: base_hmc.hpp:184
const XML_Char XML_Content * model
Definition: expat.h:151
Integrator< Hamiltonian< Model, BaseRNG > > integrator_
Definition: base_hmc.hpp:183