1 #ifndef STAN_SERVICES_INIT_INITIALIZE_STATE_HPP 2 #define STAN_SERVICES_INIT_INITIALIZE_STATE_HPP 4 #include <boost/lexical_cast.hpp> 5 #include <boost/random/additive_combine.hpp> 6 #include <boost/random/uniform_real_distribution.hpp> 7 #include <boost/random/variate_generator.hpp> 15 #include <stan/math/prim/mat.hpp> 44 void rm_indices_from_name(std::string& name) {
45 size_t x = name.find_first_of(
".");
46 if (std::string::npos == x)
return;
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));
76 template <
class Model>
77 bool are_all_pars_initialized(
const Model& model,
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;
87 template <
class Model>
88 bool validate_unconstrained_initialization(Eigen::VectorXd& cont_params,
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);
96 std::stringstream msg;
97 msg << param_names[n] <<
" initialized to invalid value (" 98 << cont_params[n] <<
")";
100 throw std::invalid_argument(msg.str());
116 template <
class Model>
122 validate_unconstrained_initialization(cont_params, model);
123 }
catch (
const std::exception& e) {
128 double init_log_prob;
129 Eigen::VectorXd init_grad = Eigen::VectorXd::Zero(model.num_params_r());
133 }
catch (
const std::exception& e) {
135 writer(
"Rejecting initial value:");
136 writer(
" Error evaluating the log probability " 137 "at the initial value.");
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.");
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 " 154 writer(
" Stan can't start sampling from this initial value.");
172 template <
class Model>
176 cont_params.setZero();
194 template <
class Model,
class RNG>
196 Eigen::VectorXd& cont_params,
200 int num_init_tries = -1;
202 boost::random::uniform_real_distribution<double>
203 init_range_distribution(-R, R);
205 boost::variate_generator
206 <RNG&, boost::random::uniform_real_distribution<double> >
207 init_rng(base_rng, init_range_distribution);
210 static int MAX_INIT_TRIES = 100;
212 for (num_init_tries = 1; num_init_tries <= MAX_INIT_TRIES;
214 for (
int i = 0; i < cont_params.size(); ++i)
215 cont_params(i) = init_rng();
220 if (num_init_tries > MAX_INIT_TRIES) {
221 std::stringstream R_ss, MAX_INIT_TRIES_ss;
223 MAX_INIT_TRIES_ss << MAX_INIT_TRIES;
227 writer(
"Initialization between (-" + R_ss.str() +
", " + R_ss.str()
229 + MAX_INIT_TRIES_ss.str() +
" attempts. ");
230 writer(
" Try specifying initial values," 231 " reducing ranges of constrained values," 232 " or reparameterizing the model.");
257 template <
class ContextFactory,
class Model,
class RNG>
260 Eigen::VectorXd& cont_params,
264 ContextFactory& context_factory) {
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);
272 static int MAX_INIT_TRIES = 100;
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;
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];
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);
294 cont_vecs_constrained,
297 model.transform_inits(cvc, cont_params, 0);
302 if (num_init_tries > MAX_INIT_TRIES) {
303 std::stringstream R_ss, MAX_INIT_TRIES_ss;
305 MAX_INIT_TRIES_ss << MAX_INIT_TRIES;
309 writer(
"Initialization between (-" + R_ss.str() +
", " + R_ss.str()
311 + MAX_INIT_TRIES_ss.str() +
" attempts. ");
312 writer(
" Try specifying initial values," 313 " reducing ranges of constrained values," 314 " or reparameterizing the model.");
317 }
catch(
const std::exception& e) {
318 writer(
"Initialization partially from source failed.");
346 template <
class ContextFactory,
class Model,
class RNG>
348 Eigen::VectorXd& cont_params,
353 ContextFactory& context_factory,
354 bool enable_random_init =
false,
357 typename ContextFactory::var_context_t context
358 = context_factory(source);
360 if (enable_random_init && !are_all_pars_initialized(model, context)) {
369 model.transform_inits(context, cont_params, 0);
370 }
catch(
const std::exception& e) {
371 writer(
"Initialization from source failed.");
389 val = boost::lexical_cast<
double>(s);
390 }
catch (
const boost::bad_lexical_cast& e) {
391 val = std::numeric_limits<double>::quiet_NaN();
417 template <
class ContextFactory,
class Model,
class RNG>
419 Eigen::VectorXd& cont_params,
423 ContextFactory& context_factory,
424 bool enable_random_init =
false,
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)
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...
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.
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...