Stan  2.10.0
probability, sampling & optimization
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 
7 #include <boost/math/special_functions/fpclassify.hpp>
8 #include <boost/random/variate_generator.hpp>
9 #include <boost/random/uniform_01.hpp>
10 #include <cmath>
11 #include <limits>
12 #include <stdexcept>
13 #include <string>
14 #include <vector>
15 
16 namespace stan {
17  namespace mcmc {
18 
19  template <class Model,
20  template<class, class> class Hamiltonian,
21  template<class> class Integrator,
22  class BaseRNG>
23  class base_hmc : public base_mcmc {
24  public:
25  base_hmc(const Model &model, BaseRNG& rng)
26  : base_mcmc(),
27  z_(model.num_params_r()),
28  integrator_(),
29  hamiltonian_(model),
30  rand_int_(rng),
32  nom_epsilon_(0.1),
34  epsilon_jitter_(0.0) {}
35 
36  void
38  std::stringstream nominal_stepsize;
39  nominal_stepsize << "Step size = " << get_nominal_stepsize();
40  writer(nominal_stepsize.str());
41  z_.write_metric(writer);
42  }
43 
44  void get_sampler_diagnostic_names(std::vector<std::string>& model_names,
45  std::vector<std::string>& names) {
46  z_.get_param_names(model_names, names);
47  }
48 
49  void get_sampler_diagnostics(std::vector<double>& values) {
50  z_.get_params(values);
51  }
52 
53  void seed(const Eigen::VectorXd& q) {
54  z_.q = q;
55  }
56 
57  void
60  this->hamiltonian_.init(this->z_, info_writer, error_writer);
61  }
62 
63  void
66  ps_point z_init(this->z_);
67 
68  // Skip initialization for extreme step sizes
69  if (this->nom_epsilon_ == 0 || this->nom_epsilon_ > 1e7)
70  return;
71 
72  this->hamiltonian_.sample_p(this->z_, this->rand_int_);
73  this->hamiltonian_.init(this->z_, info_writer, error_writer);
74 
75  // Guaranteed to be finite if randomly initialized
76  double H0 = this->hamiltonian_.H(this->z_);
77 
78  this->integrator_.evolve(this->z_, this->hamiltonian_,
79  this->nom_epsilon_,
80  info_writer, error_writer);
81 
82  double h = this->hamiltonian_.H(this->z_);
83  if (boost::math::isnan(h))
84  h = std::numeric_limits<double>::infinity();
85 
86  double delta_H = H0 - h;
87 
88  int direction = delta_H > std::log(0.8) ? 1 : -1;
89 
90  while (1) {
91  this->z_.ps_point::operator=(z_init);
92 
93  this->hamiltonian_.sample_p(this->z_, this->rand_int_);
94  this->hamiltonian_.init(this->z_, info_writer, error_writer);
95 
96  double H0 = this->hamiltonian_.H(this->z_);
97 
98  this->integrator_.evolve(this->z_, this->hamiltonian_,
99  this->nom_epsilon_,
100  info_writer, error_writer);
101 
102  double h = this->hamiltonian_.H(this->z_);
103  if (boost::math::isnan(h))
104  h = std::numeric_limits<double>::infinity();
105 
106  double delta_H = H0 - h;
107 
108  if ((direction == 1) && !(delta_H > std::log(0.8)))
109  break;
110  else if ((direction == -1) && !(delta_H < std::log(0.8)))
111  break;
112  else
113  this->nom_epsilon_
114  = direction == 1
115  ? 2.0 * this->nom_epsilon_
116  : 0.5 * this->nom_epsilon_;
117 
118  if (this->nom_epsilon_ > 1e7)
119  throw std::runtime_error("Posterior is improper. "
120  "Please check your model.");
121  if (this->nom_epsilon_ == 0)
122  throw std::runtime_error("No acceptably small step size could "
123  "be found. Perhaps the posterior is "
124  "not continuous?");
125  }
126 
127  this->z_.ps_point::operator=(z_init);
128  }
129 
130  typename Hamiltonian<Model, BaseRNG>::PointType& z() {
131  return z_;
132  }
133 
134  virtual void set_nominal_stepsize(double e) {
135  if (e > 0)
136  nom_epsilon_ = e;
137  }
138 
140  return this->nom_epsilon_;
141  }
142 
144  return this->epsilon_;
145  }
146 
147  virtual void set_stepsize_jitter(double j) {
148  if (j > 0 && j < 1)
149  epsilon_jitter_ = j;
150  }
151 
153  return this->epsilon_jitter_;
154  }
155 
157  this->epsilon_ = this->nom_epsilon_;
158  if (this->epsilon_jitter_)
159  this->epsilon_ *= 1.0
160  + this->epsilon_jitter_ * (2.0 * this->rand_uniform_() - 1.0);
161  }
162 
163  protected:
164  typename Hamiltonian<Model, BaseRNG>::PointType z_;
165  Integrator<Hamiltonian<Model, BaseRNG> > integrator_;
166  Hamiltonian<Model, BaseRNG> hamiltonian_;
167 
168  BaseRNG& rand_int_;
169 
170  // Uniform(0, 1) RNG
171  boost::uniform_01<BaseRNG&> rand_uniform_;
172 
173  double nom_epsilon_;
174  double epsilon_;
176  };
177 
178  } // mcmc
179 } // stan
180 #endif
Hamiltonian< Model, BaseRNG >::PointType z_
Definition: base_hmc.hpp:164
void get_sampler_diagnostics(std::vector< double > &values)
Definition: base_hmc.hpp:49
void init_hamiltonian(interface_callbacks::writer::base_writer &info_writer, interface_callbacks::writer::base_writer &error_writer)
Definition: base_hmc.hpp:58
Probability, optimization and sampling library.
virtual void set_nominal_stepsize(double e)
Definition: base_hmc.hpp:134
Point in a generic phase space.
Definition: ps_point.hpp:17
void get_sampler_diagnostic_names(std::vector< std::string > &model_names, std::vector< std::string > &names)
Definition: base_hmc.hpp:44
virtual void set_stepsize_jitter(double j)
Definition: base_hmc.hpp:147
base_writer is an abstract base class defining the interface for Stan writer callbacks.
Definition: base_writer.hpp:20
double get_current_stepsize()
Definition: base_hmc.hpp:143
void seed(const Eigen::VectorXd &q)
Definition: base_hmc.hpp:53
base_hmc(const Model &model, BaseRNG &rng)
Definition: base_hmc.hpp:25
double get_stepsize_jitter()
Definition: base_hmc.hpp:152
double get_nominal_stepsize()
Definition: base_hmc.hpp:139
boost::uniform_01< BaseRNG & > rand_uniform_
Definition: base_hmc.hpp:171
Hamiltonian< Model, BaseRNG >::PointType & z()
Definition: base_hmc.hpp:130
Hamiltonian< Model, BaseRNG > hamiltonian_
Definition: base_hmc.hpp:166
Integrator< Hamiltonian< Model, BaseRNG > > integrator_
Definition: base_hmc.hpp:165
void init_stepsize(interface_callbacks::writer::base_writer &info_writer, interface_callbacks::writer::base_writer &error_writer)
Definition: base_hmc.hpp:64
void write_sampler_state(interface_callbacks::writer::base_writer &writer)
Definition: base_hmc.hpp:37

     [ Stan Home Page ] © 2011–2016, Stan Development Team.