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.
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<PinvarUpdater> 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 << "PinvarUpdater being created" << std::endl;
clear();
_name = "Proportion of Invariable Sites";
assert(asrv);
_asrv = asrv;
}
inline PinvarUpdater::~PinvarUpdater() {
//std::cout << "PinvarUpdater being destroyed" << std::endl;
_asrv.reset();
}
inline void PinvarUpdater::clear() {
Updater::clear();
_prev_point = 0.0;
_asrv = nullptr;
reset();
}
inline double PinvarUpdater::getCurrentPoint() const {
return *(_asrv->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 > 0.0 && curr_point < 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->setPinvar(_prev_point);
}
inline void PinvarUpdater::proposeNewState() {
if (_lambda > 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->uniform();
if (p < 0.0)
p = -p;
else if (p > 1.0)
p = 1.0 - (p - 1.0);
_asrv->setPinvar(p);
_log_hastings_ratio = 0.0; //symmetric proposal
// This proposal invalidates all transition matrices and partials
_tree_manipulator->selectAllPartials();
_tree_manipulator->selectAllTMatrices();
}
}
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<OmegaUpdater> 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 << "OmegaUpdater being created" << std::endl;
clear();
_name = "Omega";
assert(q);
_q = q;
}
inline OmegaUpdater::~OmegaUpdater() {
//std::cout << "OmegaUpdater being destroyed" << std::endl;
_q.reset();
}
inline void OmegaUpdater::clear() {
Updater::clear();
_prev_point = 0.0;
_q = nullptr;
reset();
}
inline double OmegaUpdater::getCurrentPoint() const {
return *(_q->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 > 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->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->uniform() - 0.5));
_q->setOmega(m*_prev_point);
// Calculate log of Hastings ratio
_log_hastings_ratio = log(m);
// This proposal invalidates all transition matrices and partials
_tree_manipulator->selectAllPartials();
_tree_manipulator->selectAllTMatrices();
}
}