# HG changeset patch -- Bitbucket.org # Project ScytheMCMC # URL http://bitbucket.org/tristanz/scythemcmc/overview # User Tristan Zajonc # Date 1247512964 14400 # Node ID 9844e0828091518ae7a5be065db6dcd9615ada97 # Parent fd1ece9248241496a90334028751ee2ece71df9d Updates. Switched to Lecuyer and enable seed saving. --- a/ScytheBugs.txt +++ b/ScytheBugs.txt @@ -12,8 +12,23 @@ Bugs: Compiling the help using the latest doxygen seems to correct the problem, although macros don't get documented really. +- Docs for scythe::min refer to max. + +- What is linesearch? Docs are not very clear. + +- Reading in from file is fragile. outsheet in stata on OSX does not + use the proper newline characters. newlines across dos/unix/osx differ. + +- col_interchange doesn't sound like it is what it is in help file. It swaps + rows rather than reordering. + +- qnorm1 has return type "double double sd double" in docs. + Missing / Potential Features: +- rtnorm is parameterized using the variance whereas rnorm is parameterized + with sd. Possibility of errors. + - Fast versions of functions proportional to log density. In most MCMC applications sum(log(scythe::dnorm)), for instance, takes up a significant amount of time. However we usually just need the log of something @@ -64,4 +79,3 @@ Missing / Potential Features: etc. Then add a -DHAS_MKL flag. See: http://www.intel.com/software/products/mkl/docs/WebHelp/mkl.htm - --- a/src/mcmc.h +++ b/src/mcmc.h @@ -88,16 +88,23 @@ #include #include #include +#include #include #include #include #include #include +#include +#include #include #include #include #include "Simple/SimpleOpt.h" #include "Simple/SimpleIni.h" +// #include +// #include + +#include "scythe_extensions.h" /// Scythe MCMC version. #define SCYTHE_MCMC_VERSION "0.1" @@ -118,9 +125,9 @@ typedef pair ipair; double dInf = std::numeric_limits::infinity(); double iInf = std::numeric_limits::infinity(); -/// Global Mersenne-Twister random number generator. -mersenne myrng; - +/// Global L'Ecuyer random number generator. +lecuyer myrng; +//mersenne myrng; /*! \brief Basic MCMC options. * @@ -137,7 +144,7 @@ struct MCMCOptions { /// Chains. int chains; /// Seed for random number generator. - int random_seed; + unsigned long random_seed [6]; /// Config file (.ini style). string config_file; /// Out file to save parameter chains. @@ -189,6 +196,18 @@ inline string to_csv (const T& t) { return s; } +/*! \brief Convert any streaming type to string + * + */ +template +inline string to_string (const T& t) { + std::stringstream ss; + ss << t; + string s = ss.str(); + return s; +} + + /*! \brief Shows command line usage. Diplays: @@ -202,11 +221,11 @@ inline string to_csv (const T& t) { --chains=arg number of chains --sample-size=arg retained sample size --burnin=arg burn in period - --thin=arg thinning iterval - --random-seed=arg random number seed + --thin=arg thinning iterval (1 = no thinning) + --random-seed=arg random number seed (0 uses current timestamp) - Example: ./base_model --config-file=base_model.cfg --out-file=out.txt --chains=1 --sample-size=1000 - --burnin=0 --thin=1 --random-seed=12345 + Example: ./command --config-file=command.ini --out-file=out.txt --chains=1 --sample-size=1000 + --burnin=0 --thin=1 --random-seed=0 */ void ShowUsage() { @@ -217,11 +236,11 @@ void ShowUsage() { << " --chains=arg \t\t number of chains" << endl << " --sample-size=arg \t\t retained sample size" << endl << " --burnin=arg \t\t burn in period" << endl - << " --thin=arg \t\t\t thinning iterval (0 = no thinning)" << endl - << " --random-seed=arg \t\t random number seed" << endl << endl; + << " --thin=arg \t\t\t thinning iterval (1 = no thinning)" << endl + << " --random-seed=arg \t\t random number seed (0 uses current timestamp, 1 attempts to load from lecuyer.seed file)" << endl << endl; - cout << " Example: ./base_model --config-file=base_model.cfg --out-file=out.txt --chains=1 --sample-size=1000" << endl - << "\t --burnin=0 --thin=1 --random-seed=12345" << endl << endl; + cout << " Example: ./command --config-file=command.ini --out-file=out.txt --chains=1 --sample-size=1000" << endl + << "\t --burnin=0 --thin=1 --random-seed=0" << endl << endl; } /*! \brief Parses command line options. @@ -285,7 +304,36 @@ void DefaultStartUp(int argc, char* argv break; case OPT_RANDOM_SEED: setflags++; - mcmc_options.random_seed = atoi(args.OptionArg()); + unsigned long seed = atoi(args.OptionArg()); + cout << "TEST" << seed< seedmat("lecuyer.seed"); + for(int i = 0; i < 6; ++i) { + mcmc_options.random_seed[i] = seedmat(i); + } + } + else { + cout << "No seed exists..." << endl; + seed = time(NULL); + for(int i = 0; i < 6; ++i) { + mcmc_options.random_seed[i] = seed; + } + } + } + else { + for(int i = 0; i < 6; ++i) { + mcmc_options.random_seed[i] = seed; + } + } break; } } else { @@ -309,8 +357,8 @@ void DefaultStartUp(int argc, char* argv cout << "Burn in cannot be negative.\n"; exit(1); } - if (mcmc_options.thin < 0) { - cout << "Thinning interval cannot be negative.\n"; + if (mcmc_options.thin < 1) { + cout << "Thinning interval must be positive.\n"; exit(1); } if (!FileExists(mcmc_options.config_file)) { @@ -319,8 +367,9 @@ void DefaultStartUp(int argc, char* argv } // Seed random number generator. - myrng.initialize(mcmc_options.random_seed); - + //myrng.initialize(mcmc_options.random_seed); + myrng.SetSeed(mcmc_options.random_seed); + // Print mcmc options. cout << "Arguments:" << endl; cout << " Config file: " << mcmc_options.config_file << endl @@ -328,7 +377,13 @@ void DefaultStartUp(int argc, char* argv << " Desired sample size: " << mcmc_options.sample_size << endl << " Burn in period: " << mcmc_options.burnin << endl << " Thinning interval: " << mcmc_options.thin << endl - << " Random number seed: " << mcmc_options.random_seed << endl << endl; + << " Random number seed: " + << mcmc_options.random_seed[0] << " " + << mcmc_options.random_seed[1] << " " + << mcmc_options.random_seed[2] << " " + << mcmc_options.random_seed[3] << " " + << mcmc_options.random_seed[4] << " " + << mcmc_options.random_seed[5] << endl << endl; } @@ -455,7 +510,7 @@ public: /// Save a new value of the parameter. /// \param new_value New value to save. - virtual void Save(double new_value) {}; + virtual void Save(ParameterStorageType new_value) {}; /// Parameter is tracked / saved. /// \return True if parameter is tracked. @@ -505,7 +560,7 @@ public: * * Executes the main step function, such as taking a MH step or deterministic step. */ - virtual void DoStep() {} + virtual void DoStep() {} /*! \brief Set starting value. * * Calls the parameter's Parameter::Starting function and then saves the result using Parameter::Save. @@ -562,7 +617,7 @@ public: * Draws from the parameter's Parameter::RandomPosterior function and then saves the result using Parameter::Save. */ void DoStep() { - parameter_.Save(parameter_.RandomPosterior()); + parameter_.Save(parameter_.RandomPosterior()); } void Start() { @@ -614,12 +669,9 @@ public: /*! \brief Set function step starting value. * - * \note Because the starting value uses the normal function calculation, it is important to add - * parameters in an order that makes the function feasible. StartingValue functions are called in the order - * they are added. Any implemention of a StartingValue method for a FunctionStep will not be used. */ void Start() { - DoStep(); + parameter_.StartingValue(); } string ParameterLabel() { @@ -815,7 +867,7 @@ public: * * Calculates the value of Parameter::Function and then saves the result using Parameter::Save. */ - void DoStep() { + void DoStep() { // Draw a new parameter double new_value = proposal_.Draw(parameter_.Value()); double old_value = parameter_.Value(); @@ -878,14 +930,16 @@ public: * \param lower Lower bound on the support of the distribution (default = -Infinity) * \param upper Upper bound on the support of the distribution (default = Infinity) */ - SliceStep(ParameterType parameter, double w = 1.0, double lower = -dInf, double upper = dInf) : - parameter_(parameter), w_(w), lower_(lower), upper_(upper) {} + SliceStep(ParameterType parameter, double w, double lower, double upper) : + parameter_(parameter), w_(w), lower_(lower), upper_(upper) { + step_count_ = 0; + } /*! \brief Take a slice sampling step for the parameter. * */ - void DoStep() { - // Find the log density at the initial point. + void DoStep() { + // Find the log density at the initial point. double x0 = parameter_.Value(); double gx0 = parameter_.LogDensity(x0); @@ -945,7 +999,10 @@ public: L = x1; } } - + + // Adapt + //step_count_++; + //w_ = (step_count_/(step_count_+1)) * w_ + (1/(step_count_+1))*(R - L); // Save the point sampled parameter_.Save(x1); } @@ -972,6 +1029,7 @@ private: double w_; double lower_; double upper_; + int step_count_; }; /*! \brief MCMC sampler. @@ -1015,23 +1073,25 @@ public: * \param number_of_iterations Number of iterations to run sampler. * \param progress Display a progress bar. */ - void Iterate(int number_of_iterations, bool progress = false) { - if (progress) { - boost::progress_display show_progress(number_of_iterations); - for (int iter = 0; iter < number_of_iterations; ++iter) { - for (int i = 0; i < steps_.size(); ++i) { - steps_[i].DoStep(); - } - ++show_progress; - } - } else { - for (int iter = 0; iter < number_of_iterations; ++iter) { - for (int i = 0; i < steps_.size(); ++i) { - steps_[i].DoStep(); - } - } - } - } + void Iterate(int number_of_iterations, bool progress = false) { + if(progress) { + boost::progress_display show_progress(number_of_iterations); + for(int iter = 0; iter < number_of_iterations; ++iter) { + for(int i = 0; i < steps_.size(); ++i) { + steps_[i].DoStep(); + } + ++show_progress; + } + } + else { + for(int iter = 0; iter < number_of_iterations; ++iter) { + for(int i = 0; i < steps_.size(); ++i) { + steps_[i].DoStep(); + } + } + } + } + /*! \brief Run sampler. * * Run the MCMC sampling procedure, including burning in, sampling, thinning, displaying progress, @@ -1075,7 +1135,7 @@ public: cout << endl << "Sampling..." << endl; boost::progress_display show_progress(sample_size_); for (int i = 0; i < sample_size_; ++i) { - Iterate(thin_ + 1); + Iterate(thin_); // Save data for (int k = 0; k < tracks_.size(); ++k) { if (k == tracks_.size() - 1) { @@ -1089,7 +1149,18 @@ public: } outfile.close(); - + + // Save lecuyer seed for future use + Matrix seedmat(6, 1, false); + unsigned long seed[6]; + myrng.GetState(seed); + cout << "Saving seed: " << endl << seed[0] << " " << seed[1] << " " << seed[2] + << " " << seed[3] << " " << seed[4] << " " << seed[5] << endl; + for(int i = 0; i < 6; ++i) { + seedmat(i) = seed[i]; + } + seedmat.save("lecuyer.seed", 'o'); + cout << "Total elapsed time: " << timer.elapsed() << " seconds" << endl; } /*! \brief Number of steps in one sampler iteration. @@ -1142,3 +1213,4 @@ double fast_log_dnorm(const double x, co } +