Stan  2.10.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>
8 #include <stan/model/util.hpp>
15 #include <stan/math/prim/scal/fun/is_inf.hpp>
16 #include <stan/math/prim/scal/fun/is_nan.hpp>
17 #include <stan/math/prim/mat/fun/Eigen.hpp>
18 #include <cmath>
19 #include <iostream>
20 #include <limits>
21 #include <stdexcept>
22 #include <string>
23 #include <vector>
24 
25 namespace stan {
26  namespace services {
27  namespace init {
28  namespace {
29 
46  void rm_indices_from_name(std::string& name) {
47  size_t x = name.find_first_of(".");
48  if (std::string::npos == x) return;
49  name.erase(x);
50  }
51 
57  void rm_indices_from_name(std::vector<std::string>& names) {
58  for (size_t i = 0; i < names.size(); i++)
59  rm_indices_from_name(names[i]);
60  std::vector<std::string>::iterator it;
61  it = std::unique(names.begin(), names.end());
62  names.resize(std::distance(names.begin(), it));
63  }
64 
78  template <class Model>
79  bool are_all_pars_initialized(const Model& model,
80  const stan::io::var_context& context) {
81  std::vector<std::string> names;
82  model.constrained_param_names(names, false, false);
83  rm_indices_from_name(names);
84  for (size_t i = 0; i < names.size(); i++)
85  if (!context.contains_r(names[i])) return false;
86  return true;
87  }
88 
89  template <class Model>
90  bool validate_unconstrained_initialization(Eigen::VectorXd& cont_params,
91  Model& model) {
92  for (int n = 0; n < cont_params.size(); n++) {
93  if (stan::math::is_inf(cont_params[n])
94  || stan::math::is_nan(cont_params[n])) {
95  std::vector<std::string> param_names;
96  model.unconstrained_param_names(param_names);
97 
98  std::stringstream msg;
99  msg << param_names[n] << " initialized to invalid value ("
100  << cont_params[n] << ")";
101 
102  throw std::invalid_argument(msg.str());
103  }
104  }
105  return true;
106  }
107  }
108 
109  /***
110  * Set initial values to what container cont_params has.
111  *
112  * @param[in] cont_params the initialized state. This should be the
113  * right size and set to 0.
114  * @param[in,out] model the model. Side effects on model? I'm not
115  * quite sure
116  * @param[in,out] writer writer callback for messages
117  */
118  template <class Model>
119  bool initialize_state_values(Eigen::VectorXd& cont_params,
120  Model& model,
122  writer) {
123  try {
124  validate_unconstrained_initialization(cont_params, model);
125  } catch (const std::exception& e) {
126  writer(e.what());
127  return false;
128  }
129  double init_log_prob;
130  Eigen::VectorXd init_grad = Eigen::VectorXd::Zero(model.num_params_r());
131  try {
132  stan::model::gradient(model, cont_params, init_log_prob,
133  init_grad);
134  } catch (const std::exception& e) {
135  io::write_error_msg(writer, e);
136  writer();
137  writer("Rejecting initial value:");
138  writer(" Error evaluating the log probability "
139  "at the initial value.");
140  return false;
141  }
142  if (!boost::math::isfinite(init_log_prob)) {
143  writer("Rejecting initial value:");
144  writer(" Log probability evaluates to log(0), "
145  "i.e. negative infinity.");
146  writer(" Stan can't start sampling from this initial value.");
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  return false;
156  }
157  }
158  return true;
159  }
160 
161 
171  template <class Model>
172  bool initialize_state_zero(Eigen::VectorXd& cont_params,
173  Model& model,
175  cont_params.setZero();
176  return initialize_state_values(cont_params, model, writer);
177  }
178 
179 
193  template <class Model, class RNG>
194  bool initialize_state_random(const double R,
195  Eigen::VectorXd& cont_params,
196  Model& model,
197  RNG& base_rng,
199  int num_init_tries = -1;
200 
201  boost::random::uniform_real_distribution<double>
202  init_range_distribution(-R, R);
203 
204  boost::variate_generator
205  <RNG&, boost::random::uniform_real_distribution<double> >
206  init_rng(base_rng, init_range_distribution);
207 
208  // Random initializations until log_prob is finite
209  static int MAX_INIT_TRIES = 100;
210 
211  for (num_init_tries = 1; num_init_tries <= MAX_INIT_TRIES;
212  ++num_init_tries) {
213  for (int i = 0; i < cont_params.size(); ++i)
214  cont_params(i) = init_rng();
215  if (initialize_state_values(cont_params, model, writer))
216  break;
217  }
218 
219  if (num_init_tries > MAX_INIT_TRIES) {
220  std::stringstream R_ss, MAX_INIT_TRIES_ss;
221  R_ss << R;
222  MAX_INIT_TRIES_ss << MAX_INIT_TRIES;
223 
224  writer();
225  writer();
226  writer("Initialization between (-" + R_ss.str() + ", " + R_ss.str()
227  + ") failed after "
228  + MAX_INIT_TRIES_ss.str() + " attempts. ");
229  writer(" Try specifying initial values,"
230  " reducing ranges of constrained values,"
231  " or reparameterizing the model.");
232  return false;
233  }
234  return true;
235  }
236 
237 
256  template <class ContextFactory, class Model, class RNG>
257  bool initialize_state_source_and_random(const std::string& source,
258  double R,
259  Eigen::VectorXd& cont_params,
260  Model& model,
261  RNG& base_rng,
263  ContextFactory& context_factory) {
264  try {
265  boost::random::uniform_real_distribution<double>
266  init_range_distribution(-R, R);
267  boost::variate_generator
268  <RNG&, boost::random::uniform_real_distribution<double> >
269  init_rng(base_rng, init_range_distribution);
270 
271  static int MAX_INIT_TRIES = 100;
272 
273  int num_init_tries = -1;
274  std::vector<std::string> cont_names;
275  model.constrained_param_names(cont_names, false, false);
276  rm_indices_from_name(cont_names);
277  for (num_init_tries = 1; num_init_tries <= MAX_INIT_TRIES;
278  ++num_init_tries) {
279  std::vector<double> cont_vecs(cont_params.size());
280  for (int i = 0; i < cont_params.size(); ++i) {
281  cont_vecs[i] = init_rng();
282  cont_params[i] = cont_vecs[i];
283  }
284  typename ContextFactory::var_context_t context
285  = context_factory(source);
286  std::vector<double> cont_vecs_constrained;
287  std::vector<int> int_vecs;
288  model.write_array(base_rng, cont_vecs, int_vecs,
289  cont_vecs_constrained, false, false, 0);
290  std::vector<std::vector<size_t> > dims;
291  model.get_dims(dims);
292  stan::io::array_var_context random_context(cont_names,
293  cont_vecs_constrained,
294  dims);
295  stan::io::chained_var_context cvc(context, random_context);
296  model.transform_inits(cvc, cont_params, 0);
297  if (initialize_state_values(cont_params, model, writer))
298  break;
299  }
300 
301  if (num_init_tries > MAX_INIT_TRIES) {
302  std::stringstream R_ss, MAX_INIT_TRIES_ss;
303  R_ss << R;
304  MAX_INIT_TRIES_ss << MAX_INIT_TRIES;
305 
306  writer();
307  writer();
308  writer("Initialization between (-" + R_ss.str() + ", " + R_ss.str()
309  + ") failed after "
310  + MAX_INIT_TRIES_ss.str() + " attempts. ");
311  writer(" Try specifying initial values,"
312  " reducing ranges of constrained values,"
313  " or reparameterizing the model.");
314  return false;
315  }
316  } catch(const std::exception& e) {
317  writer("Initialization partially from source failed.");
318  writer(e.what());
319  return false;
320  }
321  return true;
322  }
323 
324 
345  template <class ContextFactory, class Model, class RNG>
346  bool initialize_state_source(const std::string source,
347  Eigen::VectorXd& cont_params,
348  Model& model,
349  RNG& base_rng,
351  writer,
352  ContextFactory& context_factory,
353  bool enable_random_init = false,
354  double R = 2) {
355  try {
356  typename ContextFactory::var_context_t context
357  = context_factory(source);
358 
359  if (enable_random_init && !are_all_pars_initialized(model, context)) {
361  R,
362  cont_params,
363  model,
364  base_rng,
365  writer,
366  context_factory);
367  }
368  model.transform_inits(context, cont_params, 0);
369  } catch(const std::exception& e) {
370  writer("Initialization from source failed.");
371  writer(e.what());
372  return false;
373  }
374  return initialize_state_values(cont_params, model, writer);
375  }
376 
377 
386  bool get_double_from_string(const std::string& s, double& val) {
387  try {
388  val = boost::lexical_cast<double>(s);
389  } catch (const boost::bad_lexical_cast& e) {
390  val = std::numeric_limits<double>::quiet_NaN();
391  return false;
392  }
393  return true;
394  }
395 
416  template <class ContextFactory, class Model, class RNG>
417  bool initialize_state(const std::string& init,
418  Eigen::VectorXd& cont_params,
419  Model& model,
420  RNG& base_rng,
422  ContextFactory& context_factory,
423  bool enable_random_init = false,
424  double init_r = 2) {
425  double R;
426  if (get_double_from_string(init, R)) {
427  if (R == 0) {
428  return initialize_state_zero(cont_params, model, writer);
429  } else {
430  return initialize_state_random(R, cont_params, model,
431  base_rng, writer);
432  }
433  }
434  return initialize_state_source(init, cont_params, model,
435  base_rng, writer,
436  context_factory,
437  enable_random_init,
438  init_r);
439  }
440 
441  } // init
442  } // services
443 } // stan
444 #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: util.hpp:422
virtual bool contains_r(const std::string &name) const =0
Return true if the specified variable name is defined.
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.
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.