Create a class whose objects can generate pseudorandom numbers.
MCMC involves random numbers at several levels. Proposed changes to the tree topology, edge lengths, or other substitution model parameters require random numbers, as does the decision to accept or reject a proposed change. The Lot class will generate the pseudorandom numbers that will be used to make these decisions.
The main member functions we will use are:
setSeed
Providing an integer to this function completely determines the sequence of pseudorandom numbers that will be issued by this object, allowing for exact replication of an analysis.
uniform and logUniform
These functions provide random deviates from a Uniform(0,1) distribution, with the latter function returning the logarithm of a Uniform(0,1) random deviate.
normal
Provides a standard normal, i.e. Normal(0,1), random deviate.
gamma
Provides a random deviate from a Gamma distribution with the specified shape and scale (mean equals shape*scale in this parameterization)
randint
Provides a random draw from a discrete uniform distribution such that the returned value is greater than or equal to low and less than or equal to high.
Create a new file lot.hpp and populate it with the following class declaration.
#pragma once
#include <ctime>
#include <memory>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/random/normal_distribution.hpp>
#include <boost/random/gamma_distribution.hpp>
#include <boost/random/variate_generator.hpp>
namespace strom {
class Lot {
public:
Lot();
~Lot();
void setSeed(unsigned seed);
double uniform();
int randint(int low, int high);
double normal();
double gamma(double shape, double scale);
double logUniform();
typedef std::shared_ptr<Lot> SharedPtr;
private:
typedef boost::variate_generator<boost::mt19937 &, boost::random::uniform_01<> > uniform_variate_generator_t;
typedef boost::variate_generator<boost::mt19937 &, boost::random::normal_distribution<> > normal_variate_generator_t;
typedef boost::variate_generator<boost::mt19937 &, boost::random::gamma_distribution<> > gamma_variate_generator_t;
typedef boost::variate_generator<boost::mt19937 &, boost::random::uniform_int_distribution<> > uniform_int_generator_t;
unsigned _seed;
boost::mt19937 _generator;
std::shared_ptr<uniform_variate_generator_t> _uniform_variate_generator;
std::shared_ptr<normal_variate_generator_t> _normal_variate_generator;
std::shared_ptr<gamma_variate_generator_t> _gamma_variate_generator;
std::shared_ptr<uniform_int_generator_t> _uniform_int_generator;
double _gamma_shape;
int _low;
int _high;
};
// member function bodies go here
}
All of the heavy lifting in our Lot class will be done by the Boost random library. While I will briefly describe the member functions, I refer you to the Boost random documentation for further details.
Boost random offers a choice of several pseudorandom number generators. I have chosen to use the Mersenne Twister 19937. This is why the _generator data member has type boost::mt19937. The constructor sets the seed of the Mersenne Twister using the current time. It is always a good idea to set your own seed so that an analysis can be repeated exactly if needed, so a setSeed member function is also provided (see below).
The constructor also creates 4 function objects (variate generators) that can translate the pseudorandom numbers output by the Mersenne Twister into random deviates from particular probability distributions. The 4 specialized functors are:
_uniform_variate_generator, which produces Uniform(0,1) random deviates,_normal_variate_generator, which produces Normal(0,1) random deviates,_gamma_variate_generator, which produces Gamma(_gamma_shape, 1) random deviates, and_uniform_int_generator, which produces discrete uniform deviates in the range _low to _high, inclusive.
All are stored as shared pointers and so must be dereferenced to be called as functions.The data member _gamma_shape is initialized to 1 and represents the shape parameter of the gamma distribution encapsulated by _gamma_variate_generator. It is possible to change the value of _gamma_shape using the gamma member function.
The data members _low, and _high store parameters of the discrete uniform distribution encapsulated by _uniform_int_generator. It is possible to change the values of _low and _high using the randint member function.
The destructor resets all 4 shared pointers so that the objects they point to can be destroyed.
inline Lot::Lot() : _seed(0), _gamma_shape(1.0), _low(0), _high(100) {
//std::cout << "Constructing a Lot" << std::endl;
_generator.seed(static_cast<unsigned int>(std::time(0)));
_uniform_variate_generator = std::shared_ptr<uniform_variate_generator_t>(new uniform_variate_generator_t(_generator, boost::random::uniform_01<>()));
_normal_variate_generator = std::shared_ptr<normal_variate_generator_t>(new normal_variate_generator_t(_generator, boost::random::normal_distribution<>()));
_gamma_variate_generator = std::shared_ptr<gamma_variate_generator_t>(new gamma_variate_generator_t(_generator, boost::random::gamma_distribution<>(_gamma_shape)));
_uniform_int_generator = std::shared_ptr<uniform_int_generator_t>(new uniform_int_generator_t(_generator, boost::random::uniform_int_distribution<>(_low, _high)));
}
inline Lot::~Lot() {
//std::cout << "Destroying a Lot" << std::endl;
_uniform_variate_generator.reset();
_normal_variate_generator.reset();
_gamma_variate_generator.reset();
_uniform_int_generator.reset();
}
This function allows you to set the seed used by the Lot object’s _generator. If the seed supplied is equal to 0, the seed is taken from the system clock.
inline void Lot::setSeed(unsigned seed) {
_seed = seed;
_generator.seed(_seed > 0 ? _seed : static_cast<unsigned int>(std::time(0)));
}
This function calls the _uniform_variate_generator functor to generate a draw from a Uniform(0,1) distribution.
inline double Lot::uniform() {
double u = (*_uniform_variate_generator)();
while (u <= 0.0)
u = (*_uniform_variate_generator)();
return u;
}
This function calls the _uniform_variate_generator functor to generate a draw from a Uniform(0,1) distribution, then returns the log of that value.
inline double Lot::logUniform() {
double u = (*_uniform_variate_generator)();
while (u <= 0.0)
u = (*_uniform_variate_generator)();
return std::log(u);
}
This function calls the _normal_variate_generator functor to generate a draw from a Normal(0,1) (i.e. Standard Normal) distribution.
inline double Lot::normal() {
return (*_normal_variate_generator)();
}
If the shape parameter equals _gamma_shape, this function calls the _gamma_variate_generator functor to generate a draw from a Gamma(shape,1) distribution, then returns that value multiplied by the scale parameter. If the shape parameter does not equal _gamma_shape, then _gamma_variate_generator is first set to point to a new gamma_variate_generator_t object with the new shape parameter and _gamma_shape is set to shape.
inline double Lot::gamma(double shape, double scale) {
assert(shape > 0.0);
assert(scale > 0.0);
if (shape != _gamma_shape) {
_gamma_shape = shape;
_gamma_variate_generator.reset(new gamma_variate_generator_t(_generator, boost::random::gamma_distribution<>(_gamma_shape,1)));
}
double deviate = (*_gamma_variate_generator)();
return scale*deviate;
}
If the low parameter equals _low and the high parameter equals _high, this function calls the _uniform_int_generator functor to generate a draw from a Discrete Uniform(_low,_high) distribution. If either the low or high parameter does not equal its data member counterpart, then _uniform_int_generator is first set to point to a new _uniform_int_generator_t object with the new low and high parameters, and _low and _high are set to low and high, respectively.
inline int Lot::randint(int low, int high) {
if (low != _low || high != _high) {
_low = low;
_high = high;
_uniform_int_generator.reset(new uniform_int_generator_t(_generator, boost::random::uniform_int_distribution<>(_low,_high)));
}
return (*_uniform_int_generator)();
}