Most of the work of updating the discrete gamma rate variance parameter has already been done by creating the Updater base class. Now all we need to do is to fill in the blanks by implementing the abstract functions in Updater.
Create a new header file named gamma_ratevar_updater.hpp and add the following class declaration code to it:
#pragma once
#include "model.hpp"
#include "updater.hpp"
#include "asrv.hpp"
namespace strom {
class GammaRateVarUpdater : public Updater {
public:
typedef std::shared_ptr<GammaRateVarUpdater> SharedPtr;
GammaRateVarUpdater(ASRV::SharedPtr asrv);
~GammaRateVarUpdater();
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
}
The constructor sets _name to be “Gamma Rate Variance” and calls the clear function. Note that _name is not listed as a data member of GammaRateVarUpdater in this class’s declaration above, so why are we allowed to modify its value here? The answer is that _name is a member of the Updater base class and thus is inherited by GammaRateVarUpdater. So, it really is a data member of GammaRateVarUpdater.
inline GammaRateVarUpdater::GammaRateVarUpdater(ASRV::SharedPtr asrv) {
//std::cout << "GammaRateVarUpdater being created" << std::endl;
clear();
_name = "Gamma Rate Variance";
assert(asrv);
_asrv = asrv;
}
inline GammaRateVarUpdater::~GammaRateVarUpdater() {
//std::cout << "GammaRateVarUpdater being destroyed" << std::endl;
_asrv.reset();
}
The clear function first calls the parent class’s version (Updater::clear()) and then returns data members _prev_point and _curr_point to their just-constructed state.
inline void GammaRateVarUpdater::clear() {
Updater::clear();
_prev_point = 0.0;
_asrv = nullptr;
reset();
}
This function returns the current value of the gamma rate variance parameter.
inline double GammaRateVarUpdater::getCurrentPoint() const {
return *(_asrv->getRateVarSharedPtr());
}
The calcLogPrior function is (you’ll remember) a pure virtual function that we must provide for this class in order to compile the program; there is no generic version of this function provided by the Updater parent class. The function body below calculates and returns the log of the Gamma prior probability density at _curr_point. The shape and scale parameters of the Gamma prior distribution are extracted from the _prior_parameters vector. This function body shows that the program produced by this tutorial will not allow the user to apply anything other than a Gamma prior distribution for this parameter, but the user will be able to set the parameters of that Gamma prior distribution.
inline double GammaRateVarUpdater::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 > 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 if (curr_point == 0.0) {
if (prior_a == 1.0) {
assert(prior_b > 0.0);
return -std::log(prior_b);
}
else if (prior_a > 1.0) {
log_prior = Updater::_log_zero;
}
else {
// prior_a < 1.0
log_prior = -Updater::_log_zero;
}
}
else
log_prior = Updater::_log_zero;
return log_prior;
}
In the event that the proposed value is rejected, this function can be called to copy _prev_point to _curr_point.
inline void GammaRateVarUpdater::revert() {
_asrv->setRateVar(_prev_point);
}
This is another pure virtual function in the parent class Updater that we must provide. The proposal algorithm used here centers a proposal window of width _lambda over _prev_point, then proposes a new value uniformly from within that proposal window. If the proposed value is less than zero, the proposed value is reflected back into the valid parameter space by simply multiplying the value by -1. While counterintuitive, this proposal is symmetric: the probability density of the proposed point given _prev_point equals the probability density of _prev_point given the proposed point. Thus, the Hastings Ratio is 1, and its log is zero.
inline void GammaRateVarUpdater::proposeNewState() {
// Save copy of _curr_point in case revert is necessary.
_prev_point = getCurrentPoint();
// Propose new value using window of width _lambda centered over _prev_point
double u = _lot->uniform();
double new_point = (_prev_point - _lambda/2.0) + u*_lambda;
assert(new_point != 0.0);
if (new_point < 0.0)
new_point *= -1.0;
_asrv->setRateVar(new_point);
// Calculate log of Hastings ratio
_log_hastings_ratio = 0.0; // symmetric proposal
// This proposal invalidates all transition matrices and partials
_tree_manipulator->selectAllPartials();
_tree_manipulator->selectAllTMatrices();
}