13.1 The Lot Class

(Win version)

< 13.0 | 13.1 | 13.2 >

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 &lt;ctime&gt;
#include &lt;memory&gt;
#include &lt;boost/random/mersenne_twister.hpp&gt;
#include &lt;boost/random/uniform_real.hpp&gt;
#include &lt;boost/random/normal_distribution.hpp&gt;
#include &lt;boost/random/gamma_distribution.hpp&gt;
#include &lt;boost/random/variate_generator.hpp&gt;

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&lt;Lot&gt;    SharedPtr;

        private:
        
            typedef boost::variate_generator&lt;boost::mt19937 &, boost::random::uniform_01&lt;&gt; &gt;                uniform_variate_generator_t;
            typedef boost::variate_generator&lt;boost::mt19937 &, boost::random::normal_distribution&lt;&gt; &gt;       normal_variate_generator_t;
            typedef boost::variate_generator&lt;boost::mt19937 &, boost::random::gamma_distribution&lt;&gt; &gt;        gamma_variate_generator_t;
            typedef boost::variate_generator&lt;boost::mt19937 &, boost::random::uniform_int_distribution&lt;&gt; &gt;  uniform_int_generator_t;

            unsigned                                        _seed;
            boost::mt19937                                  _generator;
            std::shared_ptr&lt;uniform_variate_generator_t&gt;    _uniform_variate_generator;
            std::shared_ptr&lt;normal_variate_generator_t&gt;     _normal_variate_generator;
            std::shared_ptr&lt;gamma_variate_generator_t&gt;      _gamma_variate_generator;
            std::shared_ptr&lt;uniform_int_generator_t&gt;        _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.

The constructor and destructor

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:

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 &lt;&lt; "Constructing a Lot" &lt;&lt; std::endl;
        _generator.seed(static_cast&lt;unsigned int&gt;(std::time(0)));
        _uniform_variate_generator = std::shared_ptr&lt;uniform_variate_generator_t&gt;(new uniform_variate_generator_t(_generator, boost::random::uniform_01&lt;&gt;()));
        _normal_variate_generator = std::shared_ptr&lt;normal_variate_generator_t&gt;(new normal_variate_generator_t(_generator, boost::random::normal_distribution&lt;&gt;()));
        _gamma_variate_generator = std::shared_ptr&lt;gamma_variate_generator_t&gt;(new gamma_variate_generator_t(_generator, boost::random::gamma_distribution&lt;&gt;(_gamma_shape)));
        _uniform_int_generator = std::shared_ptr&lt;uniform_int_generator_t&gt;(new uniform_int_generator_t(_generator, boost::random::uniform_int_distribution&lt;&gt;(_low, _high)));
    }
        
    inline Lot::~Lot() {
        //std::cout &lt;&lt; "Destroying a Lot" &lt;&lt; std::endl;
        _uniform_variate_generator.reset();
        _normal_variate_generator.reset();
        _gamma_variate_generator.reset();
        _uniform_int_generator.reset();
    }   

The setSeed member function

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 &gt; 0 ? _seed : static_cast&lt;unsigned int&gt;(std::time(0)));
    }   

The uniform member function

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 &lt;= 0.0)
            u = (*_uniform_variate_generator)();
        return u;
    }   

The logUniform member function

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 &lt;= 0.0)
            u = (*_uniform_variate_generator)();
        return std::log(u);
    }   

The normal member function

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)();
    }   

The gamma member function

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 &gt; 0.0);
        assert(scale &gt; 0.0);
        if (shape != _gamma_shape) {
            _gamma_shape = shape;
            _gamma_variate_generator.reset(new gamma_variate_generator_t(_generator, boost::random::gamma_distribution&lt;&gt;(_gamma_shape,1)));
        }
        double deviate = (*_gamma_variate_generator)();
        return scale*deviate;
    }   

The randint member function

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&lt;&gt;(_low,_high)));
        }
        return (*_uniform_int_generator)();
    }   

< 13.0 | 13.1 | 13.2 >