Stan  1.3
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
command.hpp
Go to the documentation of this file.
1 #ifndef __STAN__GM__COMMAND_HPP__
2 #define __STAN__GM__COMMAND_HPP__
3 
4 #include <fstream>
5 #include <boost/date_time/posix_time/posix_time_types.hpp>
6 #include <boost/math/special_functions/fpclassify.hpp>
7 #include <boost/random/additive_combine.hpp> // L'Ecuyer RNG
8 #include <boost/random/uniform_real_distribution.hpp>
9 #include <stan/version.hpp>
10 #include <stan/io/cmd_line.hpp>
11 #include <stan/io/dump.hpp>
12 #include <stan/math/functions/dot.hpp>
13 #include <stan/mcmc/adaptive_hmc.hpp>
14 #include <stan/mcmc/nuts.hpp>
15 #include <stan/mcmc/nuts_diag.hpp>
16 #include <stan/mcmc/nuts_nondiag.hpp>
17 #include <stan/mcmc/nuts_massgiven.hpp>
18 #include <stan/optimization/newton.hpp>
19 #include <stan/optimization/nesterov_gradient.hpp>
20 
21 namespace stan {
22 
23  namespace gm {
24 
25 
26  void print_nuts_help(std::string cmd) {
27  using stan::io::print_help_option;
28 
29  std::cout << std::endl;
30  std::cout << "Compiled Stan Graphical Model Command" << std::endl;
31  std::cout << std::endl;
32 
33  std::cout << "USAGE: " << cmd << " [options]" << std::endl;
34  std::cout << std::endl;
35 
36  std::cout << "OPTIONS:" << std::endl;
37  std::cout << std::endl;
38 
39  print_help_option(&std::cout,
40  "help","",
41  "Display this information");
42 
43  print_help_option(&std::cout,
44  "data","file",
45  "Read data from specified dump-format file",
46  "required if model declares data");
47 
48  print_help_option(&std::cout,
49  "init","file",
50  "Use initial values from specified file or zero values if <file>=0",
51  "default is random initialization");
52 
53  print_help_option(&std::cout,
54  "samples","file",
55  "File into which samples are written",
56  "default = samples.csv");
57 
58  print_help_option(&std::cout,
59  "append_samples","",
60  "Append samples to existing file if it exists",
61  "does not write header in append mode");
62 
63  print_help_option(&std::cout,
64  "seed","int",
65  "Random number generation seed",
66  "default = randomly generated from time");
67 
68  print_help_option(&std::cout,
69  "chain_id","int",
70  "Markov chain identifier",
71  "default = 1");
72 
73  print_help_option(&std::cout,
74  "iter","+int",
75  "Total number of iterations, including warmup",
76  "default = 2000");
77 
78  print_help_option(&std::cout,
79  "warmup","+int",
80  "Discard the specified number of initial samples",
81  "default = iter / 2");
82 
83  print_help_option(&std::cout,
84  "thin","+int",
85  "Period between saved samples after warm up",
86  "default = max(1, floor(iter - warmup) / 1000)");
87 
88  print_help_option(&std::cout,
89  "refresh","int",
90  "Period between samples updating progress report print (0 for no printing)",
91  "default = max(1,iter/200))");
92 
93  print_help_option(&std::cout,
94  "leapfrog_steps","int",
95  "Number of leapfrog steps; -1 for no-U-turn adaptation",
96  "default = -1");
97 
98  print_help_option(&std::cout,
99  "max_treedepth","int",
100  "Limit NUTS leapfrog steps to 2^max_tree_depth; -1 for no limit",
101  "default = 10");
102 
103  print_help_option(&std::cout,
104  "epsilon","float",
105  "Initial value for step size, or -1 to set automatically",
106  "default = -1");
107 
108  print_help_option(&std::cout,
109  "epsilon_pm","[0,1]",
110  "Sample epsilon +/- epsilon * epsilon_pm",
111  "default = 0.0");
112 
113  print_help_option(&std::cout,
114  "equal_step_sizes","",
115  "Use same step size for every parameter with NUTS",
116  "default is to estimate varying step sizes during warmup");
117 
118  print_help_option(&std::cout,
119  "delta","[0,1]",
120  "Accuracy target for step-size adaptation (higher means smaller step sizes)",
121  "default = 0.5");
122 
123  print_help_option(&std::cout,
124  "gamma","+float",
125  "Gamma parameter for dual averaging step-size adaptation",
126  "default = 0.05");
127 
128  print_help_option(&std::cout,
129  "save_warmup","",
130  "Save the warmup samples");
131 
132  print_help_option(&std::cout,
133  "test_grad","",
134  "Test gradient calculations using finite differences");
135 
136  print_help_option(&std::cout,
137  "point_estimate","",
138  "Fit point estimate of hidden parameters by maximizing log joint probability using Nesterov's accelerated gradient method");
139 
140  print_help_option(&std::cout,
141  "point_estimate_newton","",
142  "Fit point estimate of hidden parameters by maximizing log joint probability using Newton's method");
143 
144  print_help_option(&std::cout,
145  "nondiag_mass","",
146  "Use a nondiagonal matrix to do the sampling");
147 
148  print_help_option(&std::cout,
149  "cov_matrix","file",
150  "Preset an estimated covariance matrix");
151 
152  std::cout << std::endl;
153  }
154 
155  bool do_print(int n, int refresh) {
156  return (refresh > 0)
157  && (n == 0
158  || ((n + 1) % refresh == 0) );
159  }
160 
161  template <class Sampler, class Model, class RNG>
162  void sample_from(Sampler& sampler,
163  bool epsilon_adapt,
164  int refresh,
165  int num_iterations,
166  int num_warmup,
167  int num_thin,
168  bool save_warmup,
169  std::ostream& sample_file_stream,
170  std::vector<double>& params_r,
171  std::vector<int>& params_i,
172  Model& model,
173  RNG& base_rng) {
174 
175  sampler.set_params(params_r,params_i);
176 
177  int it_print_width = std::ceil(std::log10(num_iterations));
178  std::cout << std::endl;
179 
180  if (epsilon_adapt)
181  sampler.adapt_on();
182  for (int m = 0; m < num_iterations; ++m) {
183  if (do_print(m,refresh)) {
184  std::cout << "Iteration: ";
185  std::cout << std::setw(it_print_width) << (m + 1)
186  << " / " << num_iterations;
187  std::cout << " [" << std::setw(3)
188  << static_cast<int>((100.0 * (m + 1))/num_iterations)
189  << "%] ";
190  std::cout << ((m < num_warmup) ? " (Adapting)" : " (Sampling)");
191  std::cout << std::endl;
192  std::cout.flush();
193  }
194  if (m < num_warmup) {
195  if (save_warmup && (m % num_thin) == 0) {
196  stan::mcmc::sample sample = sampler.next();
197 
198  // FIXME: use csv_writer arg to make comma optional?
199  sample_file_stream << sample.log_prob() << ',';
200  sampler.write_sampler_params(sample_file_stream);
201  sample.params_r(params_r);
202  sample.params_i(params_i);
203  model.write_csv(base_rng,params_r,params_i,sample_file_stream,&std::cout);
204  } else {
205  sampler.next(); // discard
206  }
207  } else {
208  if (epsilon_adapt && sampler.adapting()) {
209  sampler.adapt_off();
210  sampler.write_adaptation_params(sample_file_stream);
211  }
212  if (((m - num_warmup) % num_thin) != 0) {
213  sampler.next();
214  continue;
215  } else {
216  stan::mcmc::sample sample = sampler.next();
217 
218  // FIXME: use csv_writer arg to make comma optional?
219  sample_file_stream << sample.log_prob() << ',';
220  sampler.write_sampler_params(sample_file_stream);
221  sample.params_r(params_r);
222  sample.params_i(params_i);
223  model.write_csv(base_rng,params_r,params_i,sample_file_stream,&std::cout);
224  }
225  }
226  }
227  }
228 
229  void write_comment(std::ostream& o) {
230  o << "#" << std::endl;
231  }
232  template <typename M>
233  void write_comment(std::ostream& o,
234  const M& msg) {
235  o << "# " << msg << std::endl;
236  }
237  template <typename K, typename V>
238  void write_comment_property(std::ostream& o,
239  const K& key,
240  const V& val) {
241  o << "# " << key << "=" << val << std::endl;
242  }
243 
244  template <class Model>
245  int nuts_command(int argc, const char* argv[]) {
246 
247  stan::io::cmd_line command(argc,argv);
248 
249  if (command.has_flag("help")) {
250  print_nuts_help(argv[0]);
251  return 0;
252  }
253 
254  std::string data_file;
255  command.val("data",data_file);
256  std::fstream data_stream(data_file.c_str(),
257  std::fstream::in);
258  stan::io::dump data_var_context(data_stream);
259  data_stream.close();
260 
261  Model model(data_var_context, &std::cout);
262 
263  bool point_estimate = command.has_flag("point_estimate");
264  bool point_estimate_newton = command.has_flag("point_estimate_newton");
265 
266  std::string sample_file = "samples.csv";
267  command.val("samples",sample_file);
268 
269  unsigned int num_iterations = 2000U;
270  command.val("iter",num_iterations);
271 
272  unsigned int num_warmup = num_iterations / 2;
273  command.val("warmup",num_warmup);
274 
275  unsigned int calculated_thin = (num_iterations - num_warmup) / 1000U;
276  unsigned int num_thin = (calculated_thin > 1) ? calculated_thin : 1U;
277  command.val("thin",num_thin);
278 
279  bool user_supplied_thin = command.has_key("thin");
280 
281  int leapfrog_steps = -1;
282  command.val("leapfrog_steps",leapfrog_steps);
283 
284  double epsilon = -1.0;
285  command.val("epsilon",epsilon);
286 
287  int max_treedepth = 10;
288  command.val("max_treedepth",max_treedepth);
289 
290  double epsilon_pm = 0.0;
291  command.val("epsilon_pm",epsilon_pm);
292  if (epsilon_pm < 0.0 || epsilon_pm > 1.0) {
293  std::stringstream ss;
294  ss << "epsilon_pm must be between 0 and 1"
295  << "; found epsilon_pm=" << epsilon_pm;
296  throw std::invalid_argument(ss.str());
297  }
298 
299  bool epsilon_adapt = epsilon <= 0.0;
300 
301  bool equal_step_sizes = command.has_flag("equal_step_sizes");
302 
303  double delta = 0.5;
304  command.val("delta", delta);
305 
306  double gamma = 0.05;
307  command.val("gamma", gamma);
308 
309  int refresh = num_iterations / 200;
310  refresh = refresh <= 0 ? 1 : refresh; // just for default
311  command.val("refresh",refresh);
312 
313  bool nondiag_mass = command.has_flag("nondiag_mass");
314 
315  std::string cov_file = "";
316  command.val("cov_matrix", cov_file);
317 
318  unsigned int random_seed = 0;
319  if (command.has_key("seed")) {
320  bool well_formed = command.val("seed",random_seed);
321  if (!well_formed) {
322  std::string seed_val;
323  command.val("seed",seed_val);
324  std::cerr << "value for seed must be integer"
325  << "; found value=" << seed_val << std::endl;
326  return -1;
327  }
328  } else {
329  random_seed
330  = (boost::posix_time::microsec_clock::universal_time() -
331  boost::posix_time::ptime(boost::posix_time::min_date_time))
332  .total_milliseconds();
333  }
334 
335  int chain_id = 1;
336  if (command.has_key("chain_id")) {
337  bool well_formed = command.val("chain_id",chain_id);
338  if (!well_formed || chain_id < 0) {
339  std::string chain_id_val;
340  command.val("chain_id",chain_id_val);
341  std::cerr << "value for chain_id must be positive integer"
342  << "; found chain_id=" << chain_id_val
343  << std::endl;
344  return -1;
345  }
346  }
347 
348  // FASTER, but no parallel guarantees:
349  // typedef boost::mt19937 rng_t;
350  // rng_t base_rng(static_cast<unsigned int>(random_seed + chain_id - 1);
351 
352  typedef boost::ecuyer1988 rng_t;
353  rng_t base_rng(random_seed);
354  // (2**50 = 1T samples, 1000 chains)
355  static boost::uintmax_t DISCARD_STRIDE = static_cast<boost::uintmax_t>(1) << 50;
356  // DISCARD_STRIDE <<= 50;
357  base_rng.discard(DISCARD_STRIDE * (chain_id - 1));
358 
359  std::vector<int> params_i;
360  std::vector<double> params_r;
361 
362  std::string init_val;
363  // parameter initialization
364  int num_init_tries = 1; // up here for printing below
365  if (command.has_key("init")) {
366  num_init_tries = -1;
367  command.val("init",init_val);
368  if (init_val == "0") {
369  params_i = std::vector<int>(model.num_params_i(),0);
370  params_r = std::vector<double>(model.num_params_r(),0.0);
371  } else {
372  try {
373  std::fstream init_stream(init_val.c_str(),std::fstream::in);
374  if (init_stream.fail()) {
375  std::string msg("ERROR: specified init file does not exist: ");
376  msg += init_val;
377  throw std::invalid_argument(msg);
378  }
379  stan::io::dump init_var_context(init_stream);
380  init_stream.close();
381  model.transform_inits(init_var_context,params_i,params_r);
382  } catch (const std::exception& e) {
383  std::cerr << "Error during user-specified initialization:"
384  << std::endl
385  << e.what()
386  << std::endl;
387  return -5;
388  }
389  }
390  } else {
391  init_val = "random initialization"; // for I/O
392  // init_rng generates uniformly from -2 to 2
393  boost::random::uniform_real_distribution<double>
394  init_range_distribution(-2.0,2.0);
395  boost::variate_generator<rng_t&,
396  boost::random::uniform_real_distribution<double> >
397  init_rng(base_rng,init_range_distribution);
398 
399  params_i = std::vector<int>(model.num_params_i(),0);
400  params_r = std::vector<double>(model.num_params_r());
401 
402  // retry inits until get a finite log prob value
403  std::vector<double> init_grad;
404  static int MAX_INIT_TRIES = 100;
405  for (num_init_tries = 1; num_init_tries <= MAX_INIT_TRIES; ++num_init_tries) {
406  for (size_t i = 0; i < params_r.size(); ++i)
407  params_r[i] = init_rng();
408  // FIXME: allow config vs. std::cout
409  double init_log_prob;
410  try {
411  init_log_prob = model.grad_log_prob(params_r,params_i,init_grad,&std::cout);
412  } catch (std::domain_error e) {
413  stan::mcmc::write_error_msgs(&std::cout, e);
414  init_log_prob = -std::numeric_limits<double>::infinity();
415  }
416  if (!boost::math::isfinite(init_log_prob))
417  continue;
418  for (size_t i = 0; i < init_grad.size(); ++i)
419  if (!boost::math::isfinite(init_grad[i]))
420  continue;
421  break;
422  }
423  if (num_init_tries > MAX_INIT_TRIES) {
424  std::cout << std::endl << std::endl
425  << "Initialization failed after " << MAX_INIT_TRIES
426  << " attempts. "
427  << " Try specifying initial values,"
428  << " reducing ranges of constrained values,"
429  << " or reparameterizing the model."
430  << std::endl;
431  return -1;
432  }
433  }
434 
435  bool save_warmup = command.has_flag("save_warmup");
436 
437  bool append_samples = command.has_flag("append_samples");
438  std::ios_base::openmode samples_append_mode
439  = append_samples
440  ? (std::fstream::out | std::fstream::app)
441  : std::fstream::out;
442 
443  if (command.has_flag("test_grad")) {
444  std::cout << std::endl << "TEST GRADIENT MODE" << std::endl;
445  return model.test_gradients(params_r,params_i);
446  }
447 
448  if (point_estimate_newton) {
449  std::cout << "STAN OPTIMIZATION COMMAND" << std::endl;
450  if (data_file == "")
451  std::cout << "data = (specified model requires no data)" << std::endl;
452  else
453  std::cout << "data = " << data_file << std::endl;
454 
455  std::cout << "init = " << init_val << std::endl;
456  if (num_init_tries > 0)
457  std::cout << "init tries = " << num_init_tries << std::endl;
458 
459  std::cout << "output = " << sample_file << std::endl;
460  std::cout << "save_warmup = " << save_warmup<< std::endl;
461 
462  std::cout << "seed = " << random_seed
463  << " (" << (command.has_key("seed")
464  ? "user specified"
465  : "randomly generated") << ")"
466  << std::endl;
467 
468  std::fstream sample_stream(sample_file.c_str(),
469  samples_append_mode);
470 
471  write_comment(sample_stream,"Point Estimate Generated by Stan");
472  write_comment(sample_stream);
473  write_comment_property(sample_stream,"stan_version_major",stan::MAJOR_VERSION);
474  write_comment_property(sample_stream,"stan_version_minor",stan::MINOR_VERSION);
475  write_comment_property(sample_stream,"stan_version_patch",stan::PATCH_VERSION);
476  write_comment_property(sample_stream,"data",data_file);
477  write_comment_property(sample_stream,"init",init_val);
478  write_comment_property(sample_stream,"save_warmup",save_warmup);
479  write_comment_property(sample_stream,"seed",random_seed);
480  write_comment(sample_stream);
481 
482  sample_stream << "lp__,"; // log probability first
483  model.write_csv_header(sample_stream);
484 
485  std::vector<double> gradient;
486  double lp;
487  try {
488  lp = model.grad_log_prob(params_r, params_i, gradient);
489  } catch (std::domain_error e) {
490  stan::mcmc::write_error_msgs(&std::cout, e);
491  lp = -std::numeric_limits<double>::infinity();
492  }
493 
494  double lastlp = lp - 1;
495  std::cout << "initial log joint probability = " << lp << std::endl;
496  int m = 0;
497  while ((lp - lastlp) / fabs(lp) > 1e-8) {
498  lastlp = lp;
499  lp = stan::optimization::newton_step(model, params_r, params_i);
500  std::cout << "Iteration ";
501  std::cout << std::setw(2) << (m + 1) << ". ";
502  std::cout << "Log joint probability = " << std::setw(10) << lp;
503  std::cout << ". Improved by " << (lp - lastlp) << ".";
504  std::cout << std::endl;
505  std::cout.flush();
506  m++;
507  // for (size_t i = 0; i < params_r.size(); i++)
508  // fprintf(stderr, "%f ", params_r[i]);
509  // fprintf(stderr, " %f (last = %f)\n", lp, lastlp);
510  if (save_warmup) {
511  sample_stream << lp << ',';
512  model.write_csv(base_rng,params_r,params_i,sample_stream);
513  }
514  }
515 
516  sample_stream << lp << ',';
517  model.write_csv(base_rng,params_r,params_i,sample_stream);
518 
519  return 0;
520  }
521 
522  if (point_estimate) {
523  std::cout << "STAN OPTIMIZATION COMMAND" << std::endl;
524  if (data_file == "")
525  std::cout << "data = (specified model requires no data)" << std::endl;
526  else
527  std::cout << "data = " << data_file << std::endl;
528 
529  std::cout << "init = " << init_val << std::endl;
530  if (num_init_tries > 0)
531  std::cout << "init tries = " << num_init_tries << std::endl;
532 
533  std::cout << "output = " << sample_file << std::endl;
534  std::cout << "save_warmup = " << save_warmup<< std::endl;
535 
536  std::cout << "seed = " << random_seed
537  << " (" << (command.has_key("seed")
538  ? "user specified"
539  : "randomly generated") << ")"
540  << std::endl;
541 
542  std::fstream sample_stream(sample_file.c_str(),
543  samples_append_mode);
544 
545  write_comment(sample_stream,"Point Estimate Generated by Stan");
546  write_comment(sample_stream);
547  write_comment_property(sample_stream,"stan_version_major",stan::MAJOR_VERSION);
548  write_comment_property(sample_stream,"stan_version_minor",stan::MINOR_VERSION);
549  write_comment_property(sample_stream,"stan_version_patch",stan::PATCH_VERSION);
550  write_comment_property(sample_stream,"data",data_file);
551  write_comment_property(sample_stream,"init",init_val);
552  write_comment_property(sample_stream,"save_warmup",save_warmup);
553  write_comment_property(sample_stream,"seed",random_seed);
554  write_comment(sample_stream);
555 
556  sample_stream << "lp__,"; // log probability first
557  model.write_csv_header(sample_stream);
558 
559  stan::optimization::NesterovGradient ng(model, params_r, params_i,
560  -1.0,&std::cout);
561  double lp = ng.logp();
562 
563  double lastlp = lp - 1;
564  std::cout << "initial log joint probability = " << lp << std::endl;
565  int m = 0;
566  for (size_t i = 0; i < num_iterations; i++) {
567  lastlp = lp;
568  lp = ng.step();
569  ng.params_r(params_r);
570  if (do_print(i, refresh)) {
571  std::cout << "Iteration ";
572  std::cout << std::setw(2) << (m + 1) << ". ";
573  std::cout << "Log joint probability = " << std::setw(10) << lp;
574  std::cout << ". Improved by " << (lp - lastlp) << ".";
575  std::cout << std::endl;
576  std::cout.flush();
577  }
578  m++;
579  if (save_warmup) {
580  sample_stream << lp << ',';
581  model.write_csv(base_rng,params_r,params_i,sample_stream);
582  }
583  }
584 
585  sample_stream << lp << ',';
586  model.write_csv(base_rng,params_r,params_i,sample_stream);
587 
588  return 0;
589  }
590 
591  std::cout << "STAN SAMPLING COMMAND" << std::endl;
592  if (data_file == "")
593  std::cout << "data = (specified model requires no data)" << std::endl;
594  else
595  std::cout << "data = " << data_file << std::endl;
596 
597  std::cout << "init = " << init_val << std::endl;
598  if (num_init_tries > 0)
599  std::cout << "init tries = " << num_init_tries << std::endl;
600 
601  std::cout << "samples = " << sample_file << std::endl;
602  std::cout << "append_samples = " << append_samples << std::endl;
603  std::cout << "save_warmup = " << save_warmup<< std::endl;
604 
605  std::cout << "seed = " << random_seed
606  << " (" << (command.has_key("seed")
607  ? "user specified"
608  : "randomly generated") << ")"
609  << std::endl;
610  std::cout << "chain_id = " << chain_id
611  << " (" << (command.has_key("chain_id")
612  ? "user specified"
613  : "default") << ")"
614  << std::endl;
615 
616  std::cout << "iter = " << num_iterations << std::endl;
617  std::cout << "warmup = " << num_warmup << std::endl;
618  std::cout << "thin = " << num_thin
619  << (user_supplied_thin ? " (user supplied)" : " (default)")
620  << std::endl;
621 
622  std::cout << "equal_step_sizes = " << equal_step_sizes << std::endl;
623  std::cout << "nondiag_mass = " << nondiag_mass << std::endl;
624  std::cout << "leapfrog_steps = " << leapfrog_steps << std::endl;
625  std::cout << "max_treedepth = " << max_treedepth << std::endl;;
626  std::cout << "epsilon = " << epsilon << std::endl;;
627  std::cout << "epsilon_pm = " << epsilon_pm << std::endl;;
628  std::cout << "delta = " << delta << std::endl;
629  std::cout << "gamma = " << gamma << std::endl;
630 
631  std::fstream sample_stream(sample_file.c_str(),
632  samples_append_mode);
633 
634  write_comment(sample_stream,"Samples Generated by Stan");
635  write_comment(sample_stream);
636  write_comment_property(sample_stream,"stan_version_major",stan::MAJOR_VERSION);
637  write_comment_property(sample_stream,"stan_version_minor",stan::MINOR_VERSION);
638  write_comment_property(sample_stream,"stan_version_patch",stan::PATCH_VERSION);
639  write_comment_property(sample_stream,"data",data_file);
640  write_comment_property(sample_stream,"init",init_val);
641  write_comment_property(sample_stream,"append_samples",append_samples);
642  write_comment_property(sample_stream,"save_warmup",save_warmup);
643  write_comment_property(sample_stream,"seed",random_seed);
644  write_comment_property(sample_stream,"chain_id",chain_id);
645  write_comment_property(sample_stream,"iter",num_iterations);
646  write_comment_property(sample_stream,"warmup",num_warmup);
647  write_comment_property(sample_stream,"thin",num_thin);
648  write_comment_property(sample_stream,"nondiag_mass", nondiag_mass);
649  write_comment_property(sample_stream,"equal_step_sizes",equal_step_sizes);
650  write_comment_property(sample_stream,"leapfrog_steps",leapfrog_steps);
651  write_comment_property(sample_stream,"max_treedepth",max_treedepth);
652  write_comment_property(sample_stream,"epsilon",epsilon);
653  write_comment_property(sample_stream,"epsilon_pm",epsilon_pm);
654  write_comment_property(sample_stream,"delta",delta);
655  write_comment_property(sample_stream,"gamma",gamma);
656  write_comment(sample_stream);
657 
658 
659  clock_t start = clock();
660  if (nondiag_mass) {
661  stan::mcmc::nuts_nondiag<rng_t> nuts_nondiag_sampler(model,params_r,params_i,
662  max_treedepth, epsilon,
663  epsilon_pm, epsilon_adapt,
664  delta, gamma,
665  base_rng);
666 
667  // cut & paste (see below) to enable sample-specific params
668  if (!append_samples) {
669  sample_stream << "lp__,"; // log probability first
670  nuts_nondiag_sampler.write_sampler_param_names(sample_stream);
671  model.write_csv_header(sample_stream);
672  }
673  nuts_nondiag_sampler.set_error_stream(std::cout); // cout intended
674  nuts_nondiag_sampler.set_output_stream(std::cout);
675 
676  sample_from(nuts_nondiag_sampler,epsilon_adapt,refresh,
677  num_iterations,num_warmup,num_thin,save_warmup,
678  sample_stream,params_r,params_i,
679  model,
680  base_rng);
681  }
682  else if (leapfrog_steps < 0 && !equal_step_sizes) {
683  // NUTS II (with varying step size estimation during warmup)
684  stan::mcmc::nuts_diag<rng_t> nuts2_sampler(model,params_r,params_i,
685  max_treedepth, epsilon,
686  epsilon_pm, epsilon_adapt,
687  delta, gamma,
688  base_rng);
689 
690  // cut & paste (see below) to enable sample-specific params
691  if (!append_samples) {
692  sample_stream << "lp__,"; // log probability first
693  nuts2_sampler.write_sampler_param_names(sample_stream);
694  model.write_csv_header(sample_stream);
695  }
696  nuts2_sampler.set_error_stream(std::cout); // cout intended
697  nuts2_sampler.set_output_stream(std::cout);
698 
699  sample_from(nuts2_sampler,epsilon_adapt,refresh,
700  num_iterations,num_warmup,num_thin,save_warmup,
701  sample_stream,params_r,params_i,
702  model,base_rng);
703 
704  } else if (leapfrog_steps < 0 && equal_step_sizes) {
705 
706  // NUTS I (equal step sizes)
707  stan::mcmc::nuts<rng_t> nuts_sampler(model,params_r,params_i,
708  max_treedepth, epsilon,
709  epsilon_pm, epsilon_adapt,
710  delta, gamma,
711  base_rng);
712 
713  nuts_sampler.set_error_stream(std::cout);
714  nuts_sampler.set_output_stream(std::cout); // cout intended
715  // cut & paste (see below) to enable sample-specific params
716  if (!append_samples) {
717  sample_stream << "lp__,"; // log probability first
718  nuts_sampler.write_sampler_param_names(sample_stream);
719  model.write_csv_header(sample_stream);
720  }
721 
722  sample_from(nuts_sampler,epsilon_adapt,refresh,
723  num_iterations,num_warmup,num_thin,save_warmup,
724  sample_stream,params_r,params_i,
725  model,base_rng);
726 
727  } else {
728 
729  // STANDARD HMC
730  stan::mcmc::adaptive_hmc<rng_t> hmc_sampler(model,params_r,params_i,
731  leapfrog_steps,
732  epsilon, epsilon_pm, epsilon_adapt,
733  delta, gamma,
734  base_rng);
735 
736  hmc_sampler.set_error_stream(std::cout); // intended
737  hmc_sampler.set_output_stream(std::cout);
738  // cut & paste (see above) to enable sample-specific params
739  if (!append_samples) {
740  sample_stream << "lp__,"; // log probability first
741  hmc_sampler.write_sampler_param_names(sample_stream);
742  model.write_csv_header(sample_stream);
743  }
744 
745  sample_from(hmc_sampler,epsilon_adapt,refresh,
746  num_iterations,num_warmup,num_thin,save_warmup,
747  sample_stream,params_r,params_i,
748  model,base_rng);
749  }
750  clock_t end = clock();
751  double deltaT = (double)(end - start) / CLOCKS_PER_SEC;
752  std::cout << std::endl
753  << "Elapsed Time: " << deltaT << " seconds"
754  << std::endl;
755 
756  sample_stream.close();
757  std::cout << std::endl << std::endl;
758  return 0;
759  }
760 
761  } // namespace prob
762 
763 
764 } // namespace stan
765 
766 #endif
     [ Stan Home Page ] © 2011–2013, Stan Development Team.