16.5 The PinvarUpdater and OmegaUpdater Classes

(Linux version)

< 16.4 | 16.5 | 16.6 >

These classes are so similar to GammaRateVarUpdater that I will simply provide the entire source code for each of these classes and make a few comments about the few differences that do exist.

The PinvarUpdater class

Other than trivial differences (e.g. the name of the class is necessarily different), this class differs from GammaRateVarUpdater primarily in the fact that a Beta prior is used for this parameter (because its support is limited to the unit interval) and the proposal is symmetric. The calcLogPrior function computes a Beta(a,b) prior density for the current point, where a and b are the two values supplied to it via Updater::setPriorParameters. A new value is proposed by drawing a value from a Uniform(curr_point - lambda/2, curr_point + lambda/2). If the proposed point is less than 0.0 or greater than 1.0, it is reflected back into the support interval. An assert at the beginning of proposeNewState ensures that _lambda is not large enough that multiple reflections are ever needed. Like the similar proposal used in the GammaRateVarUpdater, this proposal is symmetric and thus the Hastings ratio is 1 (0 on log scale).

Enter the following text into a new header file named pinvar_updater.hpp.

#pragma once    

#include "model.hpp"
#include "updater.hpp"
#include "asrv.hpp"

namespace strom {
    
    class PinvarUpdater : public Updater {
    
        public:
            typedef std::shared_ptr&lt;PinvarUpdater&gt; SharedPtr;

                                        PinvarUpdater(ASRV::SharedPtr asrv);
                                        ~PinvarUpdater();
        
            virtual void                clear();
            double                      getCurrentPoint() const;

            // mandatory overrides of pure virtual functions
            virtual double              calcLogPrior();
            virtual void                revert();
            virtual void                proposeNewState();
        
        private:
        
            double                      _prev_point;
            ASRV::SharedPtr             _asrv;
    };

    // member function bodies go here    
    inline PinvarUpdater::PinvarUpdater(ASRV::SharedPtr asrv) {
        //std::cout &lt;&lt; "PinvarUpdater being created" &lt;&lt; std::endl;
        clear();
        _name = "Proportion of Invariable Sites";
        assert(asrv);
        _asrv = asrv;
    }

    inline PinvarUpdater::~PinvarUpdater() {
        //std::cout &lt;&lt; "PinvarUpdater being destroyed" &lt;&lt; std::endl;
        _asrv.reset();
    }

    inline void PinvarUpdater::clear() {
        Updater::clear();
        _prev_point = 0.0;
        _asrv = nullptr;
        reset();
    }

    inline double PinvarUpdater::getCurrentPoint() const {
        return *(_asrv-&gt;getPinvarSharedPtr());
    }
    
    inline double PinvarUpdater::calcLogPrior() {
        // Assumes Beta(a,b) prior with mean a/(a+b) and variance a*b/((a + b + 1)*(a + b)^2)
        assert(_prior_parameters.size() == 2);
        double prior_a = _prior_parameters[0];
        double prior_b = _prior_parameters[1];
        
        double log_prior = 0.0;
        double curr_point = getCurrentPoint();
        if (curr_point &gt; 0.0 && curr_point &lt; 1.0) {
            log_prior += (prior_a - 1.0)*std::log(curr_point);
            log_prior += (prior_b - 1.0)*std::log(1.0 - curr_point);
            log_prior += std::lgamma(prior_a + prior_b);
            log_prior -= std::lgamma(prior_a);
            log_prior -= std::lgamma(prior_b);
        }
        else
            log_prior = Updater::_log_zero;
        return log_prior;
    }

    inline void PinvarUpdater::revert() {
        _asrv-&gt;setPinvar(_prev_point);
    }

    inline void PinvarUpdater::proposeNewState() {
        if (_lambda &gt; 1.0)
            _lambda = 1.0;
        
        // Save copy of _curr_point in case revert is necessary.
        _prev_point = getCurrentPoint();
        
        // Propose new value using uniform window of width _lambda centered over _prev_point
        double p = (_prev_point - _lambda/2.0) + _lambda*_lot-&gt;uniform();
        if (p &lt; 0.0)
            p = -p;
        else if (p &gt; 1.0)
            p = 1.0 - (p - 1.0);
        _asrv-&gt;setPinvar(p);
        
        _log_hastings_ratio = 0.0;  //symmetric proposal
        
        // This proposal invalidates all transition matrices and partials
        _tree_manipulator-&gt;selectAllPartials();
        _tree_manipulator-&gt;selectAllTMatrices();
    }

}   

The OmegaUpdater class

This class is essentially identical to GammaRateVarUpdater except for the fact that a shared pointer to a QMatrix object is passed into the constructor and stored instead of a shared ponter to an ASRV object. Another difference is that a multiplicative proposal (rather than a centered window proposal) is used to propose new parameter values in OmegaUpdater. The proposal algorithm used here generates a proposed value close to the old value (_prev_point) by multiplying _prev_point by a factor that is either slightly larger than or slightly smaller than 1.0. What is meant by “slightly” depends on _lambda, the boldness factor, with larger values of _lambda yielding bigger distances, on average, between _prev_point and the proposed point. This proposal is not symmetric: the probability density of proposing the new point given _prev_point is not the same as the probability density of proposing _prev_point given the new point. This proposal asymmetry means that we have to put the (log of the) ratio of these two probabilities into the _log_hastings_ratio data member so that the asymmetry of the move can be countered by a modified acceptance probability. That is, if proposals tend to be to the right twice as often as to the left, we accept half as many to the right to compensate for the proposal bias.

Enter the following text into a new header file named omega_updater.hpp.

#pragma once    

#include "model.hpp"
#include "updater.hpp"
#include "qmatrix.hpp"

namespace strom {
    
    class OmegaUpdater : public Updater {
    
        public:
            typedef std::shared_ptr&lt;OmegaUpdater&gt; SharedPtr;

                                        OmegaUpdater(QMatrix::SharedPtr q);
                                        ~OmegaUpdater();
        
            virtual void                clear();
            double                      getCurrentPoint() const;

            // mandatory overrides of pure virtual functions
            virtual double              calcLogPrior();
            virtual void                revert();
            virtual void                proposeNewState();
        
        private:
        
            double                      _prev_point;
            QMatrix::SharedPtr          _q;
    };

    // member function bodies go here    
    inline OmegaUpdater::OmegaUpdater(QMatrix::SharedPtr q) {
        //std::cout &lt;&lt; "OmegaUpdater being created" &lt;&lt; std::endl;
        clear();
        _name = "Omega";
        assert(q);
        _q = q;
    }

    inline OmegaUpdater::~OmegaUpdater() { 
        //std::cout &lt;&lt; "OmegaUpdater being destroyed" &lt;&lt; std::endl;
        _q.reset();
    }

    inline void OmegaUpdater::clear() {
        Updater::clear();
        _prev_point = 0.0;
        _q = nullptr;
        reset();
    }

    inline double OmegaUpdater::getCurrentPoint() const {
        return *(_q-&gt;getOmegaSharedPtr());
    } 
    
    inline double OmegaUpdater::calcLogPrior() {
        // Assumes Gamma(a,b) prior with mean a*b and variance a*b^2
        assert(_prior_parameters.size() == 2);
        double prior_a = _prior_parameters[0];
        double prior_b = _prior_parameters[1];
        
        double log_prior = 0.0;
        double curr_point = getCurrentPoint();
        if (curr_point &gt; 0.0) {
            log_prior += (prior_a - 1.0)*std::log(curr_point);
            log_prior -= curr_point/prior_b;
            log_prior -= prior_a*std::log(prior_b);
            log_prior -= std::lgamma(prior_a);
        }
        else
            log_prior = Updater::_log_zero;
        return log_prior;
    } 

    inline void OmegaUpdater::revert() {
        _q-&gt;setOmega(_prev_point);
    }

    inline void OmegaUpdater::proposeNewState() {
        // Save copy of _curr_point in case revert is necessary.
        _prev_point = getCurrentPoint();
        
        // Propose new value using multiplier with boldness _lambda
        double m = exp(_lambda*(_lot-&gt;uniform() - 0.5));
        _q-&gt;setOmega(m*_prev_point);
        
        // Calculate log of Hastings ratio
        _log_hastings_ratio = log(m);
        
        // This proposal invalidates all transition matrices and partials
        _tree_manipulator-&gt;selectAllPartials();
        _tree_manipulator-&gt;selectAllTMatrices();
    } 

}   

< 16.4 | 16.5 | 16.6 >