Stan  2.14.0
probability, sampling & optimization
initialize_state.hpp
Go to the documentation of this file.
1 #ifndef STAN_SERVICES_INIT_INITIALIZE_STATE_HPP
2 #define STAN_SERVICES_INIT_INITIALIZE_STATE_HPP
3 
4 #include <boost/lexical_cast.hpp>
5 #include <boost/random/additive_combine.hpp> // L'Ecuyer RNG
6 #include <boost/random/uniform_real_distribution.hpp>
7 #include <boost/random/variate_generator.hpp>
15 #include <stan/math/prim/mat.hpp>
16 #include <cmath>
17 #include <iostream>
18 #include <limits>
19 #include <stdexcept>
20 #include <string>
21 #include <vector>
22 
23 namespace stan {
24  namespace services {
25  namespace init {
26  namespace {
27 
44  void rm_indices_from_name(std::string& name) {
45  size_t x = name.find_first_of(".");
46  if (std::string::npos == x) return;
47  name.erase(x);
48  }
49 
55  void rm_indices_from_name(std::vector<std::string>& names) {
56  for (size_t i = 0; i < names.size(); i++)
57  rm_indices_from_name(names[i]);
58  std::vector<std::string>::iterator it;
59  it = std::unique(names.begin(), names.end());
60  names.resize(std::distance(names.begin(), it));
61  }
62 
76  template <class Model>
77  bool are_all_pars_initialized(const Model& model,
78  const stan::io::var_context& context) {
79  std::vector<std::string> names;
80  model.constrained_param_names(names, false, false);
81  rm_indices_from_name(names);
82  for (size_t i = 0; i < names.size(); i++)
83  if (!context.contains_r(names[i])) return false;
84  return true;
85  }
86 
87  template <class Model>
88  bool validate_unconstrained_initialization(Eigen::VectorXd& cont_params,
89  Model& model) {
90  for (int n = 0; n < cont_params.size(); n++) {
91  if (stan::math::is_inf(cont_params[n])
92  || stan::math::is_nan(cont_params[n])) {
93  std::vector<std::string> param_names;
94  model.unconstrained_param_names(param_names);
95 
96  std::stringstream msg;
97  msg << param_names[n] << " initialized to invalid value ("
98  << cont_params[n] << ")";
99 
100  throw std::invalid_argument(msg.str());
101  }
102  }
103  return true;
104  }
105  }
106 
107  /***
108  * Set initial values to what container cont_params has.
109  *
110  * @param[in] cont_params the initialized state. This should be the
111  * right size and set to 0.
112  * @param[in,out] model the model. Side effects on model? I'm not
113  * quite sure
114  * @param[in,out] writer writer callback for messages
115  */
116  template <class Model>
117  bool initialize_state_values(Eigen::VectorXd& cont_params,
118  Model& model,
120  writer) {
121  try {
122  validate_unconstrained_initialization(cont_params, model);
123  } catch (const std::exception& e) {
124  writer(e.what());
125  writer();
126  return false;
127  }
128  double init_log_prob;
129  Eigen::VectorXd init_grad = Eigen::VectorXd::Zero(model.num_params_r());
130  try {
131  stan::model::gradient(model, cont_params, init_log_prob,
132  init_grad, writer);
133  } catch (const std::exception& e) {
134  io::write_error_msg(writer, e);
135  writer("Rejecting initial value:");
136  writer(" Error evaluating the log probability "
137  "at the initial value.");
138  writer();
139  return false;
140  }
141  if (!boost::math::isfinite(init_log_prob)) {
142  writer("Rejecting initial value:");
143  writer(" Log probability evaluates to log(0), "
144  "i.e. negative infinity.");
145  writer(" Stan can't start sampling from this initial value.");
146  writer();
147  return false;
148  }
149  for (int i = 0; i < init_grad.size(); ++i) {
150  if (!boost::math::isfinite(init_grad(i))) {
151  writer("Rejecting initial value:");
152  writer(" Gradient evaluated at the initial value "
153  "is not finite.");
154  writer(" Stan can't start sampling from this initial value.");
155  writer();
156  return false;
157  }
158  }
159  return true;
160  }
161 
162 
172  template <class Model>
173  bool initialize_state_zero(Eigen::VectorXd& cont_params,
174  Model& model,
176  cont_params.setZero();
177  return initialize_state_values(cont_params, model, writer);
178  }
179 
180 
194  template <class Model, class RNG>
195  bool initialize_state_random(const double R,
196  Eigen::VectorXd& cont_params,
197  Model& model,
198  RNG& base_rng,
200  int num_init_tries = -1;
201 
202  boost::random::uniform_real_distribution<double>
203  init_range_distribution(-R, R);
204 
205  boost::variate_generator
206  <RNG&, boost::random::uniform_real_distribution<double> >
207  init_rng(base_rng, init_range_distribution);
208 
209  // Random initializations until log_prob is finite
210  static int MAX_INIT_TRIES = 100;
211 
212  for (num_init_tries = 1; num_init_tries <= MAX_INIT_TRIES;
213  ++num_init_tries) {
214  for (int i = 0; i < cont_params.size(); ++i)
215  cont_params(i) = init_rng();
216  if (initialize_state_values(cont_params, model, writer))
217  break;
218  }
219 
220  if (num_init_tries > MAX_INIT_TRIES) {
221  std::stringstream R_ss, MAX_INIT_TRIES_ss;
222  R_ss << R;
223  MAX_INIT_TRIES_ss << MAX_INIT_TRIES;
224 
225  writer();
226  writer();
227  writer("Initialization between (-" + R_ss.str() + ", " + R_ss.str()
228  + ") failed after "
229  + MAX_INIT_TRIES_ss.str() + " attempts. ");
230  writer(" Try specifying initial values,"
231  " reducing ranges of constrained values,"
232  " or reparameterizing the model.");
233  return false;
234  }
235  return true;
236  }
237 
238 
257  template <class ContextFactory, class Model, class RNG>
258  bool initialize_state_source_and_random(const std::string& source,
259  double R,
260  Eigen::VectorXd& cont_params,
261  Model& model,
262  RNG& base_rng,
264  ContextFactory& context_factory) {
265  try {
266  boost::random::uniform_real_distribution<double>
267  init_range_distribution(-R, R);
268  boost::variate_generator
269  <RNG&, boost::random::uniform_real_distribution<double> >
270  init_rng(base_rng, init_range_distribution);
271 
272  static int MAX_INIT_TRIES = 100;
273 
274  int num_init_tries = -1;
275  std::vector<std::string> cont_names;
276  model.constrained_param_names(cont_names, false, false);
277  rm_indices_from_name(cont_names);
278  for (num_init_tries = 1; num_init_tries <= MAX_INIT_TRIES;
279  ++num_init_tries) {
280  std::vector<double> cont_vecs(cont_params.size());
281  for (int i = 0; i < cont_params.size(); ++i) {
282  cont_vecs[i] = init_rng();
283  cont_params[i] = cont_vecs[i];
284  }
285  typename ContextFactory::var_context_t context
286  = context_factory(source);
287  std::vector<double> cont_vecs_constrained;
288  std::vector<int> int_vecs;
289  model.write_array(base_rng, cont_vecs, int_vecs,
290  cont_vecs_constrained, false, false, 0);
291  std::vector<std::vector<size_t> > dims;
292  model.get_dims(dims);
293  stan::io::array_var_context random_context(cont_names,
294  cont_vecs_constrained,
295  dims);
296  stan::io::chained_var_context cvc(context, random_context);
297  model.transform_inits(cvc, cont_params, 0);
298  if (initialize_state_values(cont_params, model, writer))
299  break;
300  }
301 
302  if (num_init_tries > MAX_INIT_TRIES) {
303  std::stringstream R_ss, MAX_INIT_TRIES_ss;
304  R_ss << R;
305  MAX_INIT_TRIES_ss << MAX_INIT_TRIES;
306 
307  writer();
308  writer();
309  writer("Initialization between (-" + R_ss.str() + ", " + R_ss.str()
310  + ") failed after "
311  + MAX_INIT_TRIES_ss.str() + " attempts. ");
312  writer(" Try specifying initial values,"
313  " reducing ranges of constrained values,"
314  " or reparameterizing the model.");
315  return false;
316  }
317  } catch(const std::exception& e) {
318  writer("Initialization partially from source failed.");
319  writer(e.what());
320  return false;
321  }
322  return true;
323  }
324 
325 
346  template <class ContextFactory, class Model, class RNG>
347  bool initialize_state_source(const std::string source,
348  Eigen::VectorXd& cont_params,
349  Model& model,
350  RNG& base_rng,
352  writer,
353  ContextFactory& context_factory,
354  bool enable_random_init = false,
355  double R = 2) {
356  try {
357  typename ContextFactory::var_context_t context
358  = context_factory(source);
359 
360  if (enable_random_init && !are_all_pars_initialized(model, context)) {
362  R,
363  cont_params,
364  model,
365  base_rng,
366  writer,
367  context_factory);
368  }
369  model.transform_inits(context, cont_params, 0);
370  } catch(const std::exception& e) {
371  writer("Initialization from source failed.");
372  writer(e.what());
373  return false;
374  }
375  return initialize_state_values(cont_params, model, writer);
376  }
377 
378 
387  bool get_double_from_string(const std::string& s, double& val) {
388  try {
389  val = boost::lexical_cast<double>(s);
390  } catch (const boost::bad_lexical_cast& e) {
391  val = std::numeric_limits<double>::quiet_NaN();
392  return false;
393  }
394  return true;
395  }
396 
417  template <class ContextFactory, class Model, class RNG>
418  bool initialize_state(const std::string& init,
419  Eigen::VectorXd& cont_params,
420  Model& model,
421  RNG& base_rng,
423  ContextFactory& context_factory,
424  bool enable_random_init = false,
425  double init_r = 2) {
426  double R;
427  if (get_double_from_string(init, R)) {
428  if (R == 0) {
429  return initialize_state_zero(cont_params, model, writer);
430  } else {
431  return initialize_state_random(R, cont_params, model,
432  base_rng, writer);
433  }
434  }
435  return initialize_state_source(init, cont_params, model,
436  base_rng, writer,
437  context_factory,
438  enable_random_init,
439  init_r);
440  }
441 
442  } // init
443  } // services
444 } // stan
445 #endif
void write_error_msg(interface_callbacks::writer::base_writer &writer, const std::exception &e)
Writes a Metropolis rejection message.
bool initialize_state_values(Eigen::VectorXd &cont_params, Model &model, interface_callbacks::writer::base_writer &writer)
bool initialize_state_source(const std::string source, Eigen::VectorXd &cont_params, Model &model, RNG &base_rng, interface_callbacks::writer::base_writer &writer, ContextFactory &context_factory, bool enable_random_init=false, double R=2)
Creates the initial state using the source parameter.
Probability, optimization and sampling library.
void gradient(const M &model, const Eigen::Matrix< double, Eigen::Dynamic, 1 > &x, double &f, Eigen::Matrix< double, Eigen::Dynamic, 1 > &grad_f, std::ostream *msgs=0)
Definition: gradient.hpp:14
bool initialize_state(const std::string &init, Eigen::VectorXd &cont_params, Model &model, RNG &base_rng, interface_callbacks::writer::base_writer &writer, ContextFactory &context_factory, bool enable_random_init=false, double init_r=2)
Creates the initial state.
bool get_double_from_string(const std::string &s, double &val)
Converts string to double.
A var_reader reads array variables of integer and floating point type by name and dimension...
Definition: var_context.hpp:29
bool initialize_state_random(const double R, Eigen::VectorXd &cont_params, Model &model, RNG &base_rng, interface_callbacks::writer::base_writer &writer)
Initializes state to random uniform values within range.
base_writer is an abstract base class defining the interface for Stan writer callbacks.
Definition: base_writer.hpp:20
bool initialize_state_zero(Eigen::VectorXd &cont_params, Model &model, interface_callbacks::writer::base_writer &writer)
Sets initial state to zero.
bool initialize_state_source_and_random(const std::string &source, double R, Eigen::VectorXd &cont_params, Model &model, RNG &base_rng, interface_callbacks::writer::base_writer &writer, ContextFactory &context_factory)
Creates the initial state.
virtual bool contains_r(const std::string &name) const =0
Return true if the specified variable name is defined.
A chained_var_context object represents two objects of var_context as one.
An array_var_context object represents a named arrays with dimensions constructed from an array...

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