mcmc.hpp

Go to the documentation of this file.

#ifndef _MCMC_HPP
#define _MCMC_HPP

#include "properties.hpp"
#include "random.hpp"

namespace Arak {

  template <class MarkovChainTraits>
  class MarkovChain {

  public:

    typedef typename MarkovChainTraits::StateType StateType;
    typedef typename MarkovChainTraits::DistributionType DistributionType;
    typedef typename MarkovChainTraits::ProposalType ProposalType;
    typedef typename MarkovChainTraits::MoveType MoveType;

  protected:

    Arak::Util::Random& random;  

    DistributionType& dist;

    ProposalType& proposal;

    StateType& state;

    double logLikelihood;

    unsigned long int numSteps;
  
    unsigned long int numMoves;

    unsigned int burnIn;

    unsigned int stride;

  public:

    MarkovChain(DistributionType& dist, 
    ProposalType& proposal,
    StateType& state,
    const Arak::Util::PropertyMap& props,
    Arak::Util::Random& random = Arak::Util::default_random) 
      : random(random), 
  dist(dist), 
  proposal(proposal),
  state(state), 
  logLikelihood(dist.logLikelihood(state)),
  numSteps(0),
  numMoves(0), 
  burnIn(0), 
  stride(0) { 
      using namespace Arak::Util;
      // Process the estimation parameters.
      if (hasp(props, "arak.mcmc.estimation.burn_in")) 
  assert(parse(getp(props, "arak.mcmc.estimation.burn_in"), burnIn));
      if (hasp(props, "arak.mcmc.estimation.stride")) 
  assert(parse(getp(props, "arak.mcmc.estimation.stride"), stride));
    }

    const StateType& getState() const { return state; }

    double getLogLikelihood() const { return logLikelihood; }
    
    StateType& getState() { return state; }

    const DistributionType& getDistribution() const { 
      return dist; 
    }

    virtual bool advance(){
      numSteps++;
      // Sample the proposal distribution.
      proposal.sample(state, random, true);
      MoveType& move = proposal.move();      
      // Compute the forwards and reverse proposal probabilities.
      double pll = proposal.ll(state);
      double rpll = proposal.rll(state);
      // Calculate the current log likelihood.
      double llh = dist.logLikelihood(state);
      // Perform the move and then check if it must be undone.
      move.execute(state);
      double llh_new = dist.logLikelihood(state);
      // Calculate the Metropolis-Hastings ratio.
      double mh_log_ratio = (llh_new + rpll) - (llh + pll);
      double mh_ratio = exp(mh_log_ratio);
      // Reject the move with some probability.
      if ((mh_ratio < 1.0) && (random.uniform() > mh_ratio)) {
  move.undo(state);
  proposal.result(false);
  return false;
      }
      // The move was successful.  Update the log likelihood.
      numMoves++;
      logLikelihood = dist.logLikelihood(state);
      proposal.result(true);
      return true;
    }

    unsigned long int getNumSteps() const { return numSteps; }
    
    unsigned long int getNumSamples() const { 
      unsigned long int s = (numSteps / (stride + 1));
      return (s > burnIn) ? s - burnIn : 0;
    }
    
    bool isSample() const { 
      return (getNumSamples() > burnIn) && 
  ((numSteps % (stride + 1)) == 0);
    }

    unsigned long int getNumMoves() const { return numMoves; }

    double acceptance() const { return double(numMoves) / double(numSteps);  }

  };
  
  template <class MarkovChainTraits>
  class AnnealedMarkovChain : public MarkovChain<MarkovChainTraits> {

  public:

    typedef typename MarkovChain<MarkovChainTraits>::StateType StateType;
    typedef typename MarkovChain<MarkovChainTraits>::DistributionType DistributionType;
    typedef typename MarkovChain<MarkovChainTraits>::ProposalType ProposalType;
    typedef typename MarkovChain<MarkovChainTraits>::MoveType MoveType;

  protected:

    unsigned long int duration;

    double cold;

    double beta1;

    double beta2;

    double theta;

  public:

    AnnealedMarkovChain(DistributionType& dist, 
      ProposalType& proposal,
      StateType& state,
      const Arak::Util::PropertyMap& props,
      Arak::Util::Random& random = 
      Arak::Util::default_random) 
      : MarkovChain<MarkovChainTraits>(dist, proposal, state, props, random) { 
      using namespace Arak::Util;
      // Process the annealing parameters.
      if (hasp(props, "arak.mcmc.annealing.hot")) {
  double hot;
  assert(parse(getp(props, "arak.mcmc.annealing.hot"), hot));
  assert(parse(getp(props, "arak.mcmc.annealing.cold"), cold));
  assert(parse(getp(props, "arak.mcmc.annealing.duration"), duration));
  // Compute the annealing parameters from the supplied parameters.
  if (hot != cold) {
    theta = (ln(hot) - ln(cold)) / ln(duration);
    double power = pow(double(duration), theta);
    beta1 = (cold - hot / power) * (1.0 - 1.0 / power);
    beta2 = hot - beta1;
  } else {
    beta1 = 0.0;
    beta2 = 1.0;
    theta = 0.0;
  }
      } else {
  beta1 = 1.0;
  beta2 = 0.0;
  theta = 0.0;
      }
    }

    virtual bool advance(){
      numSteps++;
      // Sample the proposal distribution.
      proposal.sample(state, random, false);
      MoveType& move = proposal.move();      
      // Compute the temperature.
      double temp;
      if (numSteps > duration)
  temp = cold;
      else
  temp = beta1 + beta2 / pow(double(numSteps), theta);
      // Compute the forwards and reverse proposal probabilities.
      double pll = proposal.ll(state);
      double rpll = proposal.rll(state);
      // Calculate the current (annealed) log likelihood.
      double llh = dist.logLikelihood(state, temp);
      // Perform the move and then check if it must be undone.
      move.execute(state);
      double llh_new = dist.logLikelihood(state, temp);
      // Calculate the Metropolis-Hastings ratio.
      double mh_log_ratio = (llh_new + rpll) - (llh + pll);
      double mh_ratio = exp(mh_log_ratio);
      // Reject the move with some probability.
      if ((mh_ratio < 1.0) && (random.uniform() > mh_ratio)) {
  move.undo(state);
  proposal.result(false);
  return false;
      }
      // The move was successful.  Update the log likelihood.
      numMoves++;
      logLikelihood = dist.logLikelihood(state);
      proposal.result(true);
      return true;
    }

  };

  template <class MarkovChainTraits>
  class StochasticHillClimber : public MarkovChain<MarkovChainTraits> {

  public:

    typedef typename MarkovChain<MarkovChainTraits>::StateType StateType;
    typedef typename MarkovChain<MarkovChainTraits>::DistributionType DistributionType;
    typedef typename MarkovChain<MarkovChainTraits>::ProposalType ProposalType;
    typedef typename MarkovChain<MarkovChainTraits>::MoveType MoveType;

    StochasticHillClimber(DistributionType& dist, 
        ProposalType& proposal,
        StateType& state,
        const Arak::Util::PropertyMap& props,
        Arak::Util::Random& random = Arak::Util::default_random) 
      : MarkovChain<MarkovChainTraits>(dist, proposal, state, props, random) { }

    virtual bool advance(){
      numSteps++;
      // Sample the proposal distribution.
      proposal.sample(state, random, false);
      MoveType& move = proposal.move();      
      // Perform the move and then check if it must be undone.
      move.execute(state);
      double llh_new = dist.logLikelihood(state);
      // Reject the move if it reduced the log likelihood.
      if (llh_new <= logLikelihood) {
  move.undo(state);
  proposal.result(false);
  return false;
      }
      // The move was successful.  Update the log likelihood.
      numMoves++;
      logLikelihood = llh_new;
      proposal.result(true);
      return true;
    }

  };

}

#endif

---