# HG changeset patch -- Bitbucket.org # Project ScytheMCMC # URL http://bitbucket.org/tristanz/scythemcmc/overview # User Tristan Zajonc # Date 1245693686 14400 # Node ID fd1ece9248241496a90334028751ee2ece71df9d # Parent 6ceff113f5e039fa2fd0567b864f12dd2256d244 Autoformatting --- a/src/mcmc.h +++ b/src/mcmc.h @@ -15,15 +15,15 @@ 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, MCMC++, but is based - on the Scythe, which provides convenient matrix types, random number generators, - and probability distributions. Scythe MCMC is different in focus from + on the Scythe, which provides convenient matrix types, random number generators, + and probability distributions. Scythe MCMC is different in focus from MCMCPack, 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 @@ -32,7 +32,7 @@ development of samplers for many models. Even though implementing models in ScytheMCMC requires more work than - implementing comparable models in Bugs + implementing comparable models in Bugs or JAGS, 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 @@ -52,9 +52,9 @@ - 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 @@ -108,9 +108,9 @@ using namespace scythe; using namespace std; /// Scythe column-oriented double matrix typedef. -typedef Matrix matrix; +typedef Matrix matrix; /// Scythe column-oriented integer matrix typedef. -typedef Matrix imatrix; +typedef Matrix imatrix; /// Pair typedef for integers. Used to store (x, y) coordinates of a matrix location. typedef pair ipair; @@ -119,7 +119,7 @@ double dInf = std::numeric_limits::infinity(); /// Global Mersenne-Twister random number generator. -mersenne myrng; +mersenne myrng; /*! \brief Basic MCMC options. @@ -129,7 +129,7 @@ mersenne myrng; */ struct MCMCOptions { /// Sample size. (iterations - burnin)/(thin + 1). - int sample_size; + int sample_size; /// Thinning interval int thin; /// Burn in period. @@ -159,8 +159,8 @@ bool FileExists(string strFilename) { int intStat; // Attempt to get the file attributes - intStat = stat(strFilename.c_str(),&stFileInfo); - if(intStat == 0) { + intStat = stat(strFilename.c_str(), &stFileInfo); + if (intStat == 0) { // We were able to get the file attributes // so the file obviously exists. blnReturn = true; @@ -173,7 +173,7 @@ bool FileExists(string strFilename) { // more details on why stat failed. blnReturn = false; } - + return(blnReturn); } @@ -181,12 +181,11 @@ bool FileExists(string strFilename) { * */ template -inline string to_csv (const T& t) -{ +inline string to_csv (const T& t) { std::stringstream ss; ss << t; string s = ss.str(); - replace(s.begin(),s.end(),' ', ','); + replace(s.begin(), s.end(), ' ', ','); return s; } @@ -209,7 +208,7 @@ inline string to_csv (const T& t) Example: ./base_model --config-file=base_model.cfg --out-file=out.txt --chains=1 --sample-size=1000 --burnin=0 --thin=1 --random-seed=12345 -*/ +*/ void ShowUsage() { cout << " Usage:" << endl << " --help \t\t\t displays help message" << endl @@ -220,7 +219,7 @@ void ShowUsage() { << " --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; } @@ -235,7 +234,7 @@ void ShowUsage() { */ 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 }; @@ -255,71 +254,70 @@ void DefaultStartUp(int argc, char* argv int setflags = 0; while (args.Next()) { if (args.LastError() == SO_SUCCESS) { - switch(args.OptionId()) { - case OPT_HELP: - ShowUsage(); - exit(0); - break; - case OPT_CONFIG_FILE: - setflags++; - mcmc_options.config_file = args.OptionArg(); - break; - case OPT_OUT_FILE: - setflags++; - mcmc_options.out_file = args.OptionArg(); - break; - case OPT_CHAINS: - setflags++; - mcmc_options.chains = atoi(args.OptionArg()); - break; - case OPT_SAMPLE_SIZE: - setflags++; - mcmc_options.sample_size = atoi(args.OptionArg()); - break; - case OPT_BURNIN: - setflags++; - mcmc_options.burnin = atoi(args.OptionArg()); - break; - case OPT_THIN: - setflags++; - mcmc_options.thin = atoi(args.OptionArg()); - break; - case OPT_RANDOM_SEED: - setflags++; - mcmc_options.random_seed = atoi(args.OptionArg()); - break; + switch (args.OptionId()) { + case OPT_HELP: + ShowUsage(); + exit(0); + break; + case OPT_CONFIG_FILE: + setflags++; + mcmc_options.config_file = args.OptionArg(); + break; + case OPT_OUT_FILE: + setflags++; + mcmc_options.out_file = args.OptionArg(); + break; + case OPT_CHAINS: + setflags++; + mcmc_options.chains = atoi(args.OptionArg()); + break; + case OPT_SAMPLE_SIZE: + setflags++; + mcmc_options.sample_size = atoi(args.OptionArg()); + break; + case OPT_BURNIN: + setflags++; + mcmc_options.burnin = atoi(args.OptionArg()); + break; + case OPT_THIN: + setflags++; + mcmc_options.thin = atoi(args.OptionArg()); + break; + case OPT_RANDOM_SEED: + setflags++; + mcmc_options.random_seed = atoi(args.OptionArg()); + break; } - } - else { + } else { cout << "Invalid argument: " << args.OptionText() << endl; exit(1); } } - + // Check that all required flags are set. - if(setflags != 7) { + if (setflags != 7) { ShowUsage(); exit(1); } // Check argument values. - if(mcmc_options.sample_size <= 0) { + if (mcmc_options.sample_size <= 0) { cout << "Sample size must be positive.\n"; exit(1); } - if(mcmc_options.burnin < 0) { + if (mcmc_options.burnin < 0) { cout << "Burn in cannot be negative.\n"; exit(1); } - if(mcmc_options.thin < 0) { + if (mcmc_options.thin < 0) { cout << "Thinning interval cannot be negative.\n"; exit(1); } - if(!FileExists(mcmc_options.config_file)) { + if (!FileExists(mcmc_options.config_file)) { cout << "Config file does not exist.\n"; exit(1); } - + // Seed random number generator. myrng.initialize(mcmc_options.random_seed); @@ -334,7 +332,7 @@ void DefaultStartUp(int argc, char* argv } -/*! \brief Main parameter class. +/*! \brief Main parameter class. * * The model is defined as many instances of parameter subclasses. Different Step types require different methods * to be implemented in each Parameter subclass. @@ -348,7 +346,7 @@ void DefaultStartUp(int argc, char* argv * GibbsStep: * - Parameter::StartingValue * - Parameter::RandomPosterior - * + * * Metropolis-Hastings: * - Parameter::StartingValue * - Parameter::LogDensity @@ -363,25 +361,25 @@ void DefaultStartUp(int argc, char* argv * \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. - class MeanParameter : public Parameter { + class MeanParameter : public Parameter { public: /// Default constructors that call base constructors. (REQUIRED) - MeanParameter() : Parameter() {} - MeanParameter(bool track, string name) : Parameter(track, name) {} + MeanParameter() : Parameter() {} + MeanParameter(bool track, string name) : Parameter(track, name) {} /// Draw a random value from analytical posterior for GibbsStep. double RandomPosterior() { @@ -422,39 +420,49 @@ public: * \param label Label of parameter for tracking purposes. */ Parameter(bool track, string label) : track_(track), label_(label) {} - + /// Function for deterministic nodes. /// \return Function value. - virtual ParameterStorageType Function() { return 0.0; }; - + virtual ParameterStorageType Function() { + return 0.0; + }; + /// Starting value. /// \return Starting value. - virtual ParameterStorageType StartingValue() { return 0.0; }; - + virtual ParameterStorageType StartingValue() { + return 0.0; + }; + /// Log of the probability density (plus constant) /// \param value Value of parameter to evaluate density at. /// \return Log of probability density (plus constant) - virtual double LogDensity(double value) { return 0.0; } + virtual double LogDensity(double value) { + return 0.0; + } /// Return a random draw from the posterior. /// Random draw from posterior is called by GibbsStep. /// \return Random draw from posterior. - virtual ParameterStorageType RandomPosterior() { return 0.0; } - + virtual ParameterStorageType RandomPosterior() { + return 0.0; + } + /// Value of parameter. /// \return Parameter value - virtual ParameterStorageType Value() { return 0.0; }; + virtual ParameterStorageType Value() { + return 0.0; + }; /// 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() { @@ -480,7 +488,7 @@ private: * * \ingroup steps * - * The step (Gibbs, Metropolis-Hastings, Slice, etc) describes how the sampler + * The step (Gibbs, Metropolis-Hastings, Slice, etc) describes how the sampler * calculates the next value of each parameter. The key function is Step::DoStep, which is * equivalent to the Execute function in a Command Pattern. * Samplers consists of a many of steps. @@ -496,34 +504,40 @@ public: /*! \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 */ - virtual string ParameterValue() { return ""; } - + virtual string ParameterValue() { + return ""; + } + /*! \brief Return parameter label * * \return The label of the parameter associated with the step instance. * \see Parameter::Label */ - virtual string ParameterLabel() { return ""; } + virtual string ParameterLabel() { + return ""; + } /*! Return the tracking status * * \return The tracking status of the parameter associated with the step instance. * \see Parameter::Track */ - virtual bool ParameterTrack() { return true; } + virtual bool ParameterTrack() { + return true; + } }; /*! \brief Gibbs step. @@ -531,7 +545,7 @@ public: * \ingroup steps * * A univariate Gibbs step draws from the Parameter::RandomPosterior function of a Parameter and then saves the result - * using Parameter:Save. + * using Parameter:Save. */ template class GibbsStep: public Step { @@ -556,25 +570,25 @@ public: } string ParameterLabel() { - return parameter_.Label(); + return parameter_.Label(); } string ParameterValue() { - return to_csv(parameter_.Value()); + return to_csv(parameter_.Value()); } bool ParameterTrack() { - return parameter_.Track(); + return parameter_.Track(); } private: /// Parameter associated with step instance. - ParameterType parameter_; + ParameterType parameter_; }; /*! \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. @@ -609,32 +623,32 @@ public: } string ParameterLabel() { - return parameter_.Label(); + return parameter_.Label(); } string ParameterValue() { - return to_csv(parameter_.Value()); + return to_csv(parameter_.Value()); } bool ParameterTrack() { - return parameter_.Track(); + return parameter_.Track(); } private: /// Parameter associated with step instance. - ParameterType parameter_; + ParameterType parameter_; }; /*! \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 - * and the LogDensity to calculate the proposal terms in the acceptance \f$\alpha\f$. For + * and the LogDensity to calculate the proposal terms in the acceptance \f$\alpha\f$. For * symmetric proposal distributions returning 0.0 is enough. */ - + /*! \brief Normal proposal for Metropolis-Hastings. * * \ingroup proposals @@ -649,7 +663,7 @@ public: * \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. @@ -657,7 +671,7 @@ public: 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. @@ -667,14 +681,14 @@ public: 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. @@ -704,9 +718,9 @@ public: * \param idenominator Tuning parameter equal to the inverse denominator \f$1/d\f$ of a Beta distribution. */ BetaProposal(double idenominator) { - denom_ = 1/idenominator; + denom_ = 1 / idenominator; } - + /*! \brief Draw from Beta proposal * * \param starting_value Starting value (mean of Beta distribution). @@ -716,9 +730,9 @@ public: 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. */ @@ -727,19 +741,19 @@ public: 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. * - * For convenience, the log normal proposal is parameterized using the unlogged mean (starting value) + * For convenience, the log normal proposal is parameterized using the unlogged mean (starting value) * and the (positive) log standard deviation. This differs from how scythe::dlnorm is parameterized. * */ @@ -754,7 +768,7 @@ public: * \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). @@ -763,9 +777,9 @@ public: 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. */ @@ -773,7 +787,7 @@ public: double logmean = log(starting_value); return dlnorm(new_value, logmean, logsd_); } - + private: double logsd_; }; @@ -781,7 +795,7 @@ private: /*! \brief Metropolis Hastings step * * \ingroup steps - * + * * Basic univariate Metropolis-Hastings step. Requires both a Parameter and Proposal instance. */ template @@ -796,7 +810,7 @@ public: * \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. @@ -807,26 +821,26 @@ public: double old_value = parameter_.Value(); // MH accept/reject criteria double alpha = parameter_.LogDensity(new_value) - parameter_.LogDensity(old_value) - + proposal_.LogDensity(old_value, new_value) - proposal_.LogDensity(new_value, old_value); - if(myrng.runif() < min(exp(alpha), 1.0)) { + + proposal_.LogDensity(old_value, new_value) - proposal_.LogDensity(new_value, old_value); + if (myrng.runif() < min(exp(alpha), 1.0)) { parameter_.Save(new_value); } } - + void Start() { parameter_.Save(parameter_.StartingValue()); } string ParameterLabel() { - return parameter_.Label(); + return parameter_.Label(); } string ParameterValue() { - return to_csv(parameter_.Value()); + return to_csv(parameter_.Value()); } bool ParameterTrack() { - return parameter_.Track(); + return parameter_.Track(); } private: /// Parameter associated with step instance. @@ -840,7 +854,7 @@ private: * * 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. @@ -864,9 +878,9 @@ 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 = 1.0, double lower = -dInf, double upper = dInf) : + parameter_(parameter), w_(w), lower_(lower), upper_(upper) {} + /*! \brief Take a slice sampling step for the parameter. * */ @@ -874,83 +888,82 @@ public: // 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. - while(true) { - if(L <= lower_) { + while (true) { + if (L <= lower_) { break; } - if(parameter_.LogDensity(L) <= logy) { + if (parameter_.LogDensity(L) <= logy) { break; } L = L - w_; } - - while(true) { - if(R >= upper_) { + + while (true) { + if (R >= upper_) { break; } - if(parameter_.LogDensity(R) <= logy) { + if (parameter_.LogDensity(R) <= logy) { break; } R = R + w_; } - + // Shrink interval to lower and upper bounds. - if(L < lower_) { + if (L < lower_) { L = lower_; } - if(R > upper_) { + if (R > upper_) { R = upper_; } - + // Sample from the interval, shrinking it on each rejection double x1, gx1; - while(true) { + while (true) { x1 = myrng.runif() * (R - L) + L; // Sample between L and R, uniformly. gx1 = parameter_.LogDensity(x1); - if(gx1 >= logy) { + if (gx1 >= logy) { break; } - - if(x1 > x0) { + + if (x1 > x0) { R = x1; - } - else { + } else { L = x1; } } - + // Save the point sampled parameter_.Save(x1); } - + void Start() { parameter_.Save(parameter_.StartingValue()); } string ParameterLabel() { - return parameter_.Label(); + return parameter_.Label(); } string ParameterValue() { - return to_csv(parameter_.Value()); + return to_csv(parameter_.Value()); } bool ParameterTrack() { - return parameter_.Track(); + return parameter_.Track(); } private: /// Parameter associated with step instance. @@ -993,7 +1006,7 @@ public: // Add step to step stack. steps_.push_back(step); // Add index to tracking stack, if necessary tracked. - if(steps_.back().ParameterTrack()) { + if (steps_.back().ParameterTrack()) { tracks_.push_back(steps_.size() - 1); } } @@ -1003,18 +1016,17 @@ public: * \param progress Display a progress bar. */ void Iterate(int number_of_iterations, bool progress = false) { - if(progress) { + 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) { + 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) { + } else { + for (int iter = 0; iter < number_of_iterations; ++iter) { + for (int i = 0; i < steps_.size(); ++i) { steps_[i].DoStep(); } } @@ -1028,32 +1040,31 @@ public: 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()); - if(!outfile.is_open()) { + if (!outfile.is_open()) { cerr << "ERROR: Cannot open file!"; exit(1); } - - for(int k = 0; k < tracks_.size(); ++k) { - if(k == tracks_.size() - 1) { + + for (int k = 0; k < tracks_.size(); ++k) { + if (k == tracks_.size() - 1) { outfile << steps_[tracks_[k]].ParameterLabel() << endl; - } - else { + } else { outfile << steps_[tracks_[k]].ParameterLabel() << ","; } } // Setting starting value: cout << "Setting starting values..." << endl; - for(int i = 0; i < steps_.size(); ++i) { + for (int i = 0; i < steps_.size(); ++i) { steps_[i].Start(); } @@ -1063,23 +1074,22 @@ public: cout << endl << "Sampling..." << endl; boost::progress_display show_progress(sample_size_); - for(int i = 0; i < sample_size_; ++i) { + for (int i = 0; i < sample_size_; ++i) { Iterate(thin_ + 1); // Save data - for(int k = 0; k < tracks_.size(); ++k) { - if(k == tracks_.size() - 1) { + for (int k = 0; k < tracks_.size(); ++k) { + if (k == tracks_.size() - 1) { outfile << steps_[tracks_[k]].ParameterValue() << endl; - } - else { + } else { outfile << steps_[tracks_[k]].ParameterValue() << ","; } } //Show progress ++show_progress; } - + outfile.close(); - + cout << "Total elapsed time: " << timer.elapsed() << " seconds" << endl; } /*! \brief Number of steps in one sampler iteration. @@ -1111,16 +1121,15 @@ private: * Useful if === is too strict. */ bool approx_equal(double a, double b) { - return abs(a-b) <= numeric_limits::epsilon() ? true : false; + return abs(a - b) <= numeric_limits::epsilon() ? true : false; } /// Fast sum of log of dnorms template -double sum_fast_log_dnorm(const Matrix &A, const double mean, const double sd) -{ +double sum_fast_log_dnorm(const Matrix &A, const double mean, const double sd) { double result = 0; - for(int i = 0; i < A.size(); ++i) { + for (int i = 0; i < A.size(); ++i) { // Slower part here is A[i], can this be speed up using an iterator? result += fast_log_dnorm(A[i], mean, sd); } @@ -1129,7 +1138,7 @@ double sum_fast_log_dnorm(const Matrix