code needed to write custom MCMC routines.  It also provides common MCMC step types, including Gibbs, Metropolis-Hastings,

    and Slice sampling.  Users can experiment with which sampling steps provide the best results and implement their own sampling

    steps as desired.

    Scythe MCMC was motivated with the need to combine simpler samplers, such as univariate slice sampling, with more specialized

    samplers for particular problems.  Future versions of Scythe MCMC will include common samplers for nonparametric Bayesian models,

    including a Collapsed Gibbs sampler for the Dirichlet Process, and Slice samplers for the Dirichlet Process and Indian Buffet Process.

    There are many choices for writing MCMC samplers.  Scythe MCMC is similar to, and heavily inspired by,

    <a href="http://darwin.eeb.uconn.edu/mcmc++/mcmc++.html">MCMC++</a>, but is based

    <a href="http://mcmcpack.wustl.edu/wiki/index.php/Main_Page">MCMCPack</a>, which

    is also based on Scythe and interfaces with R.  Users wishing to distribute a specific model to R users

    should consider contributing to MCMCPack directly. Unlike MCMCPack, Scythe MCMC makes no attempt to interface with R

    development of samplers for many models.

    Even though implementing models in ScytheMCMC requires more work than

    or <a href="http://www-ice.iarc.fr/~martyn/software/jags/">JAGS</a>, Scythe MCMC generally

    leads to faster samplers that can exploit the particular structure of a given problem.  It also

    allows development of models that cannot be sampled effectively using either Bugs or Jags, such

       - Collapsed gibbs samplers for Dirichlet Process,

       - Slice sampler for Dirichlet Process and Indian Buffet Process.

    - MIT License.

    \section gettingstarted GETTING STARTED

    \section license MIT LICENSE

using namespace std;

/// Scythe column-oriented double matrix typedef.

/// Scythe column-oriented integer matrix typedef.

/// Pair typedef for integers.  Used to store (x, y) coordinates of a matrix location.

typedef pair<int, int> ipair;

double iInf = std::numeric_limits<int>::infinity();

/// Global Mersenne-Twister random number generator.

/*! \brief Basic MCMC options.

 */

struct MCMCOptions {

  /// Sample size. (iterations - burnin)/(thin + 1).

  /// Thinning interval

  int thin;

  /// Burn in period.

  int intStat;

  // Attempt to get the file attributes

    // We were able to get the file attributes

    // so the file obviously exists.

    blnReturn = true;

    // more details on why stat failed.

    blnReturn = false;

  }

  return(blnReturn);

}

 *

 */

template <class T>

  std::stringstream ss;

  ss << t;

  string s = ss.str();

  return s;

}

     Example: ./base_model --config-file=base_model.cfg --out-file=out.txt --chains=1 --sample-size=1000

              --burnin=0 --thin=1 --random-seed=12345

</pre>

void ShowUsage() {

  cout << "  Usage:" << endl

       << "    --help \t\t\t displays help message" << 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;

  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;

}

*/

void DefaultStartUp(int argc, char* argv[]) {

  cout << "ScytheMCMC v. " << SCYTHE_MCMC_VERSION << endl << endl;

  // Config parser is based on SimpleOpt Library.

  // See: http://code.jellycan.com/simpleopt/

  enum { OPT_SAMPLE_SIZE, OPT_THIN, OPT_BURNIN, OPT_CHAINS, OPT_RANDOM_SEED, OPT_CONFIG_FILE, OPT_OUT_FILE, OPT_HELP };

  int setflags = 0;

  while (args.Next()) {

    if (args.LastError() == SO_SUCCESS) {

      }

      cout << "Invalid argument: " << args.OptionText() << endl;

      exit(1);

    }

  }

  // Check that all required flags are set.

    ShowUsage();

    exit(1);

  }

  // Check argument values.

    cout << "Sample size must be positive.\n";

    exit(1);

  }

    cout << "Burn in cannot be negative.\n";

    exit(1);

  }

    cout << "Thinning interval cannot be negative.\n";

    exit(1);

  }

    cout << "Config file does not exist.\n";

    exit(1);

  }

  // Seed random number generator.

  myrng.initialize(mcmc_options.random_seed);

}

 *

 * The model is defined as many instances of parameter subclasses.  Different Step types require different methods

 * to be implemented in each Parameter subclass.

 * GibbsStep:

 *   - Parameter::StartingValue

 *   - Parameter::RandomPosterior

 * Metropolis-Hastings:

 *   - Parameter::StartingValue

 *   - Parameter::LogDensity

 * \section examples Examples:

 *

 * \subsection example1 Normal mean parameter with normal prior, known standard deviation:

 \code

 // Priors

 double mean_mu = 0;

 double mean_sigma = 2;

 // Parameters

 double mean = 0; // free parameter

 double sd = 1; // known

 // Data

 matrix X('X.txt');

 // Define the mean parameter, implementing methods required for GibbsStep.

 public:

   /// Default constructors that call base constructors. (REQUIRED)

  /// Draw a random value from analytical posterior for GibbsStep.

  double RandomPosterior() {

   * \param label Label of parameter for tracking purposes.

  */

  Parameter(bool track, string label) : track_(track), label_(label) {}

  /// Function for deterministic nodes.

  /// \return Function value.

  /// Starting value.

  /// \return Starting value.

  /// Log of the probability density (plus constant)

  /// \param value Value of parameter to evaluate density at.

  /// \return Log of probability density (plus constant)

  /// Return a random draw from the posterior.

  /// Random draw from posterior is called by GibbsStep.

  /// \return Random draw from posterior.

  /// Value of parameter.

  /// \return Parameter value

  /// Save a new value of the parameter.

  /// \param new_value New value to save.

  virtual void Save(double new_value) {};

  /// Parameter is tracked / saved.

  /// \return True if parameter is tracked.

  bool Track() {

    return track_;

  }

  /// String label of parameter.

  /// \return Label of parameter, for purposes of saving output.

  string Label() {

 *

 * \ingroup steps

 *

 * calculates the next value of each parameter.  The key function is Step::DoStep, which is

 * equivalent to the Execute function in a <a href="http://en.wikipedia.org/wiki/Command_pattern">Command Pattern</a>.

 * Samplers consists of a many of steps.

  /*! \brief Take step.

   *

   * Executes the main step function, such as taking a MH step or deterministic step.

  virtual void DoStep() {}

  /*! \brief Set starting value.

   *

   *  Calls the parameter's Parameter::Starting function and then saves the result using Parameter::Save.

   */

  virtual void Start() {}

  /*! \brief Return parameter value

   *

   * \return String representation of parameter value for csv output.

   * \see Parameter::Value

   */

  /*! \brief Return parameter label

   *

   * \return The label of the parameter associated with the step instance.

   * \see Parameter::Label

   */

  /*! Return the tracking status

   *

   * \return The tracking status of the parameter associated with the step instance.

   * \see Parameter::Track

   */

};

/*! \brief Gibbs step.

 * \ingroup steps

 *

 * A univariate Gibbs step draws from the Parameter::RandomPosterior function of a Parameter and then saves the result

 */

template <typename ParameterType, typename ParameterStorageType>

class GibbsStep: public Step {

  }

  string ParameterLabel() {

  }

  string ParameterValue() {

  }

  bool ParameterTrack() {

  }

private:

  /// Parameter associated with step instance.

};

/*! \brief Deterministic function step.

 *

 * \ingroup steps

 * Deterministic nodes can be helpful to track summary statistics or to reduce computations through intermediate use

 * of sufficient statistics across multiple parameters.  FunctionStep calls the parameters Parameter::Function method

 * and saves the result using Parameter::Save.

  }

  string ParameterLabel() {

  }

  string ParameterValue() {

  }

  bool ParameterTrack() {

  }

private:

  /// Parameter associated with step instance.

};

/*! \defgroup proposals Metropolis-Hastings Proposals

 * The MetropStep Metropolis-Hastings class requires both a Parameter instance

 * and a Proposal instance.  Proposal instances define a Draw and LogDensity method.

 * The constructor should set a tuning parameter.

 *

 * MetroStep::DoStep uses the Draw method to draw new proposed values for the parameter

 * symmetric proposal distributions returning 0.0 is enough.

 */

/*! \brief Normal proposal for Metropolis-Hastings.

 *

 *  \ingroup proposals

   * \param standard_deviation Tuning parameter for normal proposal.

   */

  NormalProposal(double standard_deviation) : standard_deviation_(standard_deviation) {}

  /*! \brief Draw from normal proposal

   *

   * \param starting_value Starting value.

  double Draw(double starting_value) {

    return myrng.rnorm(starting_value, standard_deviation_);

  }

  /*! \brief Log probability of proposal new value, given starting value.

   *

   * \note Needs to be implemented but can be 0.0 if proposal is symmetric, which it is.

  double LogDensity(double new_value, double starting_value) {

    return 0.0; // Symmetric, doesn't matter.

  }

private:

  double standard_deviation_;

};

/*! \brief Beta proposal for Metropolis-Hastings.

 * \ingroup proposals

 *

 * Useful proposal type for parameters with support between 0 and 1.

   * \param idenominator Tuning parameter equal to the inverse denominator \f$1/d\f$ of a Beta distribution.

   */

  BetaProposal(double idenominator) {

  }

  /*! \brief Draw from Beta proposal

   *

   * \param starting_value Starting value (mean of Beta distribution).

    double beta = alpha - denom_;

    return myrng.rbeta(alpha, beta);

  }

  /*! \brief Log probability of proposal new value, given starting value.

   * \param starting_value Starting value.

   * \param new_value Proposed new value.

   */

    double beta = alpha - denom_;

    return dbeta(new_value, alpha, beta);

  }

private:

  double denom_;

};

/*! \brief Log normal proposal for Metropolis-Hastings.

 * \ingroup proposals

 *

 * Useful proposal type for parameters with support between 0 and positive infinity.

 *

 * and the (positive) log standard deviation.  This differs from how scythe::dlnorm is parameterized.

 *

 */

   * \param logsd Tuning parameter equal to the log standard deviation.

   */

  LogNormalProposal(double logsd) : logsd_(logsd) {}

  /*! \brief Draw from log normal proposal

   *

   * \param starting_value Starting value (unlogged mean of log normal).

    double logmean = log(starting_value);

    return myrng.rlnorm(logmean, logsd_);

  }

  /*! \brief Log probability of proposal new value, given starting value.

   * \param starting_value Starting value.

   * \param new_value Proposed new value.

   */

    double logmean = log(starting_value);

    return dlnorm(new_value, logmean, logsd_);

  }

private:

  double logsd_;

};

/*! \brief Metropolis Hastings step

 *

 * \ingroup steps

 * Basic univariate Metropolis-Hastings step.  Requires both a Parameter and Proposal instance.

 */

template <typename ParameterType, typename ProposalType>

   *  \param proposal A Proposal instance that implements Proposal::LogDensity and Proposal::Draw.

   */

  MetropStep(ParameterType parameter, ProposalType proposal) : parameter_(parameter), proposal_(proposal) {}

  /*! \brief Take a function step for the parameter.

   *

   *  Calculates the value of Parameter::Function and then saves the result using Parameter::Save.

    double old_value = parameter_.Value();

    // MH accept/reject criteria

    double alpha = parameter_.LogDensity(new_value) - parameter_.LogDensity(old_value)

      parameter_.Save(new_value);

    }

  }

  void Start() {

    parameter_.Save(parameter_.StartingValue());

  }

  string ParameterLabel() {

  }

  string ParameterValue() {

  }

  bool ParameterTrack() {

  }

private:

  /// Parameter associated with step instance.

 *

 * Basic univariate slice sampling step with stepping out and shrinkage. Performs a slice sampling update from an initial

 * point to a new point that leaves invariant the distribution with the specifified log density functions.

 * The log desnity function may return -Inf for points outside the support of the distribution.

 * If a lower and/or upper bound is specified for the support, the log density function will not be called outside

 * such limits.

   *  \param lower Lower bound on the support of the distribution (default = -Infinity)

   *  \param upper Upper bound on the support of the distribution (default = Infinity)

   */

  /*! \brief Take a slice sampling step for the parameter.

   *

   */

    // Find the log density at the initial point.

    double x0 = parameter_.Value();

    double gx0 = parameter_.LogDensity(x0);

    // Determine slice level, in log terms

    double logy = gx0 - myrng.rexp(1);

    //Find the initial interval to sample from.

    double u = myrng.runif() * w_;

    double L = x0 - u;

    double R = x0 + (w_ - u);

    // Expand the interval until its ends are outside the slice, or

    // until the limit on steps is reached.

    // allow infinite steps... could hang!

    // FIXME: Add check for too many iterations and exit with error message.

        break;

      }

        break;

      }

      L = L - w_;

    }

        break;

      }

        break;

      }

      R = R + w_;

    }

    // Shrink interval to lower and upper bounds.

      L = lower_;

    }

      R = upper_;

    }

    // Sample from the interval, shrinking it on each rejection

    double x1, gx1;

      x1 = myrng.runif() * (R - L) + L;  // Sample between L and R, uniformly.

      gx1 = parameter_.LogDensity(x1);

        break;

      }

        R = x1;

        L = x1;

      }

    }

    // Save the point sampled

    parameter_.Save(x1);

  }

  void Start() {

    parameter_.Save(parameter_.StartingValue());

  }

  string ParameterLabel() {

  }

  string ParameterValue() {

  }

  bool ParameterTrack() {

  }

private:

  /// Parameter associated with step instance.

    // Add step to step stack.

    steps_.push_back(step);

    // Add index to tracking stack, if necessary tracked.

      tracks_.push_back(steps_.size() - 1);

    }

  }

   * \param progress Display a progress bar.

   */

  void Iterate(int number_of_iterations, bool progress = false) {

      boost::progress_display show_progress(number_of_iterations);

          steps_[i].DoStep();

        }

        ++show_progress;

      }

          steps_[i].DoStep();

        }

      }

  void Run() {

    // Timer

    boost::timer timer;

    // Status of sampler...

    cout << "Running sampler..." << endl;

    cout << "Number of steps added: " << NumberOfSteps() << endl;

    cout << "Number of tracked steps added: " << NumberOfTrackedSteps() << endl;

    // Opening output file

    cout << "Opening output file: " << out_file_ << endl;

    ofstream outfile(out_file_.c_str());

      cerr << "ERROR: Cannot open file!";

      exit(1);

    }

        outfile << steps_[tracks_[k]].ParameterLabel() << endl;

        outfile << steps_[tracks_[k]].ParameterLabel() << ",";

      }

    }

    // Setting starting value:

    cout << "Setting starting values..." << endl;

      steps_[i].Start();

    }

    cout << endl << "Sampling..." << endl;

    boost::progress_display show_progress(sample_size_);

      Iterate(thin_ + 1);

      // Save data

          outfile << steps_[tracks_[k]].ParameterValue() << endl;