Stan  2.10.0
probability, sampling & optimization
base_static_uniform.hpp
Go to the documentation of this file.
1 #ifndef STAN_MCMC_HMC_UNIFORM_BASE_STATIC_UNIFORM_HPP
2 #define STAN_MCMC_HMC_UNIFORM_BASE_STATIC_UNIFORM_HPP
3 
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 {
21  template <class Model,
22  template<class, class> class Hamiltonian,
23  template<class> class Integrator,
24  class BaseRNG>
26  public base_hmc<Model, Hamiltonian, Integrator, BaseRNG> {
27  public:
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 
35 
36  sample
37  transition(sample& init_sample,
40  this->sample_stepsize();
41 
42  this->seed(init_sample.cont_params());
43 
44  this->hamiltonian_.sample_p(this->z_, this->rand_int_);
45  this->hamiltonian_.init(this->z_, info_writer, error_writer);
46 
47  ps_point z_init(this->z_);
48  double H0 = this->hamiltonian_.H(this->z_);
49 
50  ps_point z_sample(this->z_);
51  double sum_prob = 1;
52  double sum_metro_prob = 1;
53 
54  boost::random::uniform_int_distribution<> uniform(0, L_ - 1);
55  int Lp = uniform(this->rand_int_);
56 
57  for (int l = 0; l < Lp; ++l) {
58  this->integrator_.evolve(this->z_, this->hamiltonian_,
59  -this->epsilon_,
60  info_writer, error_writer);
61 
62  double h = this->hamiltonian_.H(this->z_);
63  if (boost::math::isnan(h))
64  h = std::numeric_limits<double>::infinity();
65 
66  double prob = std::exp(H0 - h);
67  sum_prob += prob;
68  sum_metro_prob += prob > 1 ? 1 : prob;
69 
70  if (this->rand_uniform_() < prob / sum_prob)
71  z_sample = this->z_;
72  }
73 
74  this->z_.ps_point::operator=(z_init);
75 
76  for (int l = 0; l < L_ - 1 - Lp; ++l) {
77  this->integrator_.evolve(this->z_, this->hamiltonian_,
78  this->epsilon_,
79  info_writer, error_writer);
80 
81  double h = this->hamiltonian_.H(this->z_);
82  if (boost::math::isnan(h))
83  h = std::numeric_limits<double>::infinity();
84 
85  double prob = std::exp(H0 - h);
86  sum_prob += prob;
87  sum_metro_prob += prob > 1 ? 1 : prob;
88 
89  if (this->rand_uniform_() < prob / sum_prob)
90  z_sample = this->z_;
91  }
92 
93  double accept_prob = sum_metro_prob / static_cast<double>(L_);
94 
95  this->z_.ps_point::operator=(z_sample);
96  this->energy_ = this->hamiltonian_.H(this->z_);
97  return sample(this->z_.q,
98  - this->hamiltonian_.V(this->z_),
99  accept_prob);
100  }
101 
102  void get_sampler_param_names(std::vector<std::string>& names) {
103  names.push_back("stepsize__");
104  names.push_back("int_time__");
105  names.push_back("energy__");
106  }
107 
108  void get_sampler_params(std::vector<double>& values) {
109  values.push_back(this->epsilon_);
110  values.push_back(this->T_);
111  values.push_back(this->energy_);
112  }
113 
114  void set_nominal_stepsize_and_T(const double e, const double t) {
115  if (e > 0 && t > 0) {
116  this->nom_epsilon_ = e;
117  T_ = t;
118  update_L_();
119  }
120  }
121 
122  void set_nominal_stepsize_and_L(const double e, const int l) {
123  if (e > 0 && l > 0) {
124  this->nom_epsilon_ = e;
125  L_ = l;
126  T_ = this->nom_epsilon_ * L_; }
127  }
128 
129  void set_T(const double t) {
130  if (t > 0) {
131  T_ = t;
132  update_L_();
133  }
134  }
135 
136  void set_nominal_stepsize(const double e) {
137  if (e > 0) {
138  this->nom_epsilon_ = e;
139  update_L_();
140  }
141  }
142 
143  double get_T() {
144  return this->T_;
145  }
146 
147  int get_L() {
148  return this->L_;
149  }
150 
151  protected:
152  double T_;
153  int L_;
154  double energy_;
155 
156  void update_L_() {
157  L_ = static_cast<int>(T_ / this->nom_epsilon_);
158  L_ = L_ < 1 ? 1 : L_;
159  }
160  };
161  } // mcmc
162 } // stan
163 #endif
Hamiltonian Monte Carlo implementation that uniformly samples from trajectories with a static integra...
Hamiltonian< Model, BaseRNG >::PointType z_
Definition: base_hmc.hpp:164
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)
Definition: sample.hpp:17
Probability, optimization and sampling library.
void get_sampler_param_names(std::vector< std::string > &names)
Point in a generic phase space.
Definition: ps_point.hpp:17
void set_nominal_stepsize_and_T(const double e, const double t)
double cont_params(int k) const
Definition: sample.hpp:25
base_static_uniform(const Model &model, BaseRNG &rng)
base_writer is an abstract base class defining the interface for Stan writer callbacks.
Definition: base_writer.hpp:20
void get_sampler_params(std::vector< double > &values)
void seed(const Eigen::VectorXd &q)
Definition: base_hmc.hpp:53
void set_nominal_stepsize_and_L(const double e, const int l)
boost::uniform_01< BaseRNG & > rand_uniform_
Definition: base_hmc.hpp:171
Hamiltonian< Model, BaseRNG > hamiltonian_
Definition: base_hmc.hpp:166
Integrator< Hamiltonian< Model, BaseRNG > > integrator_
Definition: base_hmc.hpp:165
sample transition(sample &init_sample, interface_callbacks::writer::base_writer &info_writer, interface_callbacks::writer::base_writer &error_writer)

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