dualaverage.hpp

Go to the documentation of this file.

#ifndef __STAN__MCMC__DUALAVERAGE_H__
#define __STAN__MCMC__DUALAVERAGE_H__

#include <cmath>
#include <cstddef>

#include <vector>

#include <iostream>

namespace stan {

  namespace mcmc {

    class DualAverage {
    protected:
      std::vector<double> _gbar, _xbar, _x0, _lastx;
      int _k;
      double _gamma;

    public:
      DualAverage(double gamma, const std::vector<double>& x0) 
        : _gbar(x0.size(), 0), _xbar(x0.size(), 0), _x0(x0), _lastx(x0),
          _k(0), _gamma(gamma){
      }

      void update(const std::vector<double>& g, std::vector<double>& xk) {
        _k++;
        xk.resize(_gbar.size());
        double avgeta = 1.0 / (_k + 10);
        double xbar_avgeta = pow(_k, -0.75);
        double muk = 0.5 * sqrt(_k) / _gamma;
        for (size_t i = 0; i < _gbar.size(); ++i) {
          _gbar[i] = avgeta * g[i] + (1 - avgeta) * _gbar[i];
          xk[i] = _x0[i] - muk * _gbar[i];
          // fprintf(stderr, "DUALAVERAGE update %d: g = %f, gbar = %f, lastx = %f",
          //   _k, g[0], _gbar[0], _lastx[0]);
          _lastx[i] = xk[i];
          _xbar[i] = xbar_avgeta * xk[i] + (1 - xbar_avgeta) * _xbar[i];
        }
        // fprintf(stderr, ", xk = %f\n", xk[0]);
      }

      inline void setx0(const std::vector<double>& x0) {
        _x0.assign(x0.begin(), x0.end());
      }

      inline void xbar(std::vector<double>& xbar) {
        xbar.assign(_xbar.begin(), _xbar.end());
      }
      inline void gbar(std::vector<double>& gbar) {
        gbar.assign(_gbar.begin(), _gbar.end());
      }
      inline void xk(std::vector<double>& xk) {
        xk.assign(_lastx.begin(), _lastx.end());
      }
      inline void x0(std::vector<double>& x0) {
        x0.assign(_x0.begin(), _x0.end());
      }
      inline int k() { return _k; }
      inline double gamma() { return _gamma; }
    };

  }
}

#endif