Estimation C++ DLL: Adding the new Proposal class with common, adjustable seed...

Estimation C++ DLL: Adding the new Proposal class with common, adjustable seed for both normal and uniform rng and the associated changes including the removal of now obsolete RandSampler class.

mex/sources/estimation/RandSampler.cc → mex/sources/estimation/Proposal.cc

 /*
- * Copyright (C) 2010 Dynare Team
+ * Copyright (C) 2009-2010 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -18,29 +18,74 @@
  */
 
 ///////////////////////////////////////////////////////////
-//  RandSampler.cpp
-//  Implementation of the Class RandSampler
-//  Created on:      09-Sep-2010 23:25:30
+//  Proposal.cpp
+//  Implementation of the Class Proposal
+//  Created on:      15-Dec-2010 12:43:49
 ///////////////////////////////////////////////////////////
 
-#include "RandSampler.hh"
-#include "BlasBindings.hh"
+#include "Proposal.hh"
+
+Proposal::Proposal(const VectorConstView &vJscale, const MatrixConstView &covariance) :
+  len(covariance.getCols()),
+  covarianceCholeskyDecomposition(len),
+  newDraw(len),
+  uniform_rng_type(0, 1), // uniform random number generator distribution type
+  uniformVrng(base_rng, uniform_rng_type), // uniform random variate_generator
+  normal_rng_type(0, 1), // normal random number generator distribution type (mean, standard)
+  normalVrng(base_rng, normal_rng_type) // normal random variate_generator
+{
+  Matrix Jscale(len);
+  Matrix DD(len);
+  DD = covariance;
+
+  lapack::choleskyDecomp(DD, "U");
+  Jscale.setAll(0.0);
+  for (size_t i = 0; i < len; i++)
+    Jscale(i, i) = vJscale(i);
+  blas::gemm("N", "N", 1.0, DD, Jscale, 0.0, covarianceCholeskyDecomposition);
+}
+
+void
+Proposal::draw(Vector &mean, Vector &draw)
+{
+  assert(len == draw.getSize());
+  assert(len == mean.getSize());
+
+  draw = mean;
+  for (size_t i = 0; i < len; ++i)
+    newDraw(i) =  normalVrng();
+  blas::gemv("T", 1.0, covarianceCholeskyDecomposition, newDraw, 1.0, draw);
+
+}
+
+Matrix &
+Proposal::getVar()
+{
+  return covarianceCholeskyDecomposition;
+}
 
 /**
- * draw = Mean + randn(1,n) * Sigma_upper_chol;
+ * set or get if null arguments
  */
+int
+Proposal::seed()
+{
+  return curSeed;
+}
+
 void
-RandSampler::randMultiVar(Prior &distribution, Vector &draw, const Vector &mean, const Matrix &Scale, const size_t n)
+Proposal::seed(int newSeed)
+{
+  curSeed = newSeed;
+  base_rng.seed(curSeed);
+}
+
+/**
+ * currently returns uniform for MH sampler,
+ */
+double
+Proposal::selectionTestDraw()
 {
-  assert(n == draw.getSize());
-  assert(n == mean.getSize() );
-  assert(n == Scale.getRows());
-  assert(n == Scale.getCols());
-
-  draw=mean;
-  Vector drawTmp(n);
-  for (size_t i = 0; i < n; ++i)
-    drawTmp(i) =  distribution.drand();
-  blas::gemv("T", 1.0, Scale, drawTmp, 1.0, draw);
+  return uniformVrng();
 }
 

mex/sources/estimation/Proposal.hh

+/*
+ * Copyright (C) 2009-2010 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+///////////////////////////////////////////////////////////
+//  Proposal.hh
+//  Implementation of the Class Proposal
+//  Created on:      15-Dec-2010 12:43:49
+///////////////////////////////////////////////////////////
+
+#if !defined(CABFAB46_2FA5_4178_8D1C_05DFEAFB4570__INCLUDED_)
+#define CABFAB46_2FA5_4178_8D1C_05DFEAFB4570__INCLUDED_
+
+/**
+ * Proposal class will then have the common, seed initialised base generator
+ * (currently congruental rand48) and member functions such as seed, reset (to
+ * initial,default value) and have all rand generators we need that generate from
+ * that common base: single uniform and normal (either single or multivariate
+ * determined by  the size of the variance matrix)  for now.
+ *
+ * The rng interfaces will be - a single, selectionTest() currently uniform for MH
+ * sampler, draw(mean, newDraw) - single or multivariate parameters (currently
+ * from Normal(mean,Sigma))  and poss. also the seed state (set and get), the
+ * normal Chol_variance_decomp (constructed using Chol Decomp of init Variance
+ * matrix) , as core members  .This class will handle the main boost random rng
+ * intricacies. See enclosed updated diagram (if ok I will upload it)
+ *
+ */
+
+#include "Matrix.hh"
+#include "BlasBindings.hh"
+#include "LapackBindings.hh"
+#include <boost/random/linear_congruential.hpp>
+#include <boost/random/normal_distribution.hpp>
+#include <boost/random/uniform_real.hpp>
+#include <boost/random/variate_generator.hpp>
+
+// typedef for base_uniform_generator_type
+// rend48 seems better than the basic minstd_rand  but
+// one my later try ecuyer1988 or taus88 which are better but not yet supported in the Boost versionwe use
+typedef boost::rand48 base_uniform_generator_type;
+
+class Proposal
+{
+
+public:
+  Proposal();
+
+public:
+  Proposal(const VectorConstView &vJscale, const MatrixConstView &covariance);
+  virtual ~Proposal(){};
+  virtual void draw(Vector &mean, Vector &draw);
+  virtual Matrix&getVar();
+  virtual int seed();
+  virtual void seed(int seedInit);
+  virtual double selectionTestDraw();
+
+private:
+  size_t len;
+  Matrix covarianceCholeskyDecomposition;
+  /**
+   * Vector of new draws
+   *
+   */
+  Vector newDraw;
+
+  base_uniform_generator_type base_rng;
+  boost::uniform_real<> uniform_rng_type;
+  boost::variate_generator<base_uniform_generator_type &,  boost::uniform_real<> > uniformVrng;
+  boost::normal_distribution<double> normal_rng_type;
+  boost::variate_generator<base_uniform_generator_type &, boost::normal_distribution<double> > normalVrng;
+
+  int curSeed;
+
+};
+
+#endif // !defined(CABFAB46_2FA5_4178_8D1C_05DFEAFB4570__INCLUDED_)

mex/sources/estimation/RandSampler.hh

-/*
- * Copyright (C) 2010 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-///////////////////////////////////////////////////////////
-//  RandSampler.hh
-//  Implementation of the Class RandSampler
-//  Created on:      09-Sep-2010 23:25:29
-///////////////////////////////////////////////////////////
-
-#if !defined(RS5D7E5E52_2A4F_4f98_9A5C_A7FD8C278E0A__INCLUDED_)
-#define RS5D7E5E52_2A4F_4f98_9A5C_A7FD8C278E0A__INCLUDED_
-
-#include "LogPosteriorDensity.hh"
-
-class RandSampler {
-
-public:
-  virtual ~RandSampler(){};
-  virtual double compute(VectorView &mhLogPostDens, MatrixView &mhParams, Matrix &steadyState,
-                         Vector &estParams, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H,
-                         const size_t presampleStart, int &info, const size_t startDraw, size_t nMHruns, const Matrix &Jscale,
-                         LogPosteriorDensity &logPosteriorDensity, Prior &drawDistribution,
-                         EstimatedParametersDescription &epd) = 0;
-  /**
-   * draw = Mean + randn(1,n) * Sigma_upper_chol;
-   *
-   */
-  virtual void randMultiVar(Prior &distribution, Vector &draw, const Vector &mean, const Matrix &Scale, const size_t n);
-
-};
-
-#endif // !defined(5D7E5E52_2A4F_4f98_9A5C_A7FD8C278E0A__INCLUDED_)

mex/sources/estimation/RandomWalkMetropolisHastings.cc

@@ -25,8 +25,8 @@
 double
 RandomWalkMetropolisHastings::compute(VectorView &mhLogPostDens, MatrixView &mhParams, Matrix &steadyState,
                                       Vector &estParams, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H,
-                                      const size_t presampleStart, int &info, const size_t startDraw, size_t nMHruns, const Matrix &Dscale,
-                                      LogPosteriorDensity &lpd, Prior &drawDistribution, EstimatedParametersDescription &epd)
+                                      const size_t presampleStart, int &info, const size_t startDraw, size_t nMHruns, 
+                                      LogPosteriorDensity &lpd, Proposal &pDD, EstimatedParametersDescription &epd)
 {
   //streambuf *likbuf, *drawbuf *backup;
   std::ofstream urandfilestr, drawfilestr;
@@ -43,7 +43,7 @@ RandomWalkMetropolisHastings::compute(VectorView &mhLogPostDens, MatrixView &mhP
   for (size_t run = startDraw - 1; run < nMHruns; ++run)
     {
       overbound=false;
-      randMultiVar(drawDistribution, newParDraw, parDraw, Dscale, parDraw.getSize());
+      pDD.draw(parDraw, newParDraw);
       for (count=0;count<parDraw.getSize();++count)
         {
           overbound=(newParDraw(count) <  epd.estParams[count].lower_bound || newParDraw(count) > epd.estParams[count].upper_bound );
@@ -68,7 +68,7 @@ RandomWalkMetropolisHastings::compute(VectorView &mhLogPostDens, MatrixView &mhP
               newLogpost = -INFINITY;
             }
         }
-      urand=uniform.drand();
+      urand=pDD.selectionTestDraw();
       if ((newLogpost > -INFINITY) && log(urand) < newLogpost-logpost)
         {
           parDraw = newParDraw;

mex/sources/estimation/RandomWalkMetropolisHastings.hh

@@ -26,26 +26,25 @@
 #if !defined(A6BBC5E0_598E_4863_B7FF_E87320056B80__INCLUDED_)
 #define A6BBC5E0_598E_4863_B7FF_E87320056B80__INCLUDED_
 
-#include "RandSampler.hh"
+#include "LogPosteriorDensity.hh"
+#include "Proposal.hh"
 
-class RandomWalkMetropolisHastings : public RandSampler
+class RandomWalkMetropolisHastings
 {
 
 private:
-  UniformPrior uniform;
   Vector parDraw, newParDraw;
 
 public:
   RandomWalkMetropolisHastings(size_t size) :
-    uniform(0.0, 0.0, 0.0, 1.0, 0.0, 1.0),
     parDraw(size), newParDraw(size)
   {
   };
   virtual ~RandomWalkMetropolisHastings(){};
   virtual double compute(VectorView &mhLogPostDens, MatrixView &mhParams, Matrix &steadyState,
                          Vector &estParams, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H,
-                         const size_t presampleStart, int &info, const size_t startDraw, size_t nMHruns, const Matrix &Jscale,
-                         LogPosteriorDensity &logPosteriorDensity, Prior &drawDistribution,
+                         const size_t presampleStart, int &info, const size_t startDraw, size_t nMHruns, 
+                         LogPosteriorDensity &logPosteriorDensity, Proposal &proposalDrawDistribution,
                          EstimatedParametersDescription &epd);
 };
 

mex/sources/estimation/logMHMCMCposterior.cc

@@ -104,9 +104,9 @@ fillEstParamsInfo(const mxArray *estim_params_info, EstimatedParameter::pType ty
 
 int
 sampleMHMC(LogPosteriorDensity &lpd, RandomWalkMetropolisHastings &rwmh,
-           Matrix &steadyState, Vector &estParams, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H,
-           size_t presampleStart, int &info, const VectorConstView &nruns, size_t fblock, size_t nBlocks,
-           const Matrix &Jscale, Prior &drawDistribution, EstimatedParametersDescription &epd,
+           Matrix &steadyState, Vector &estParams, Vector &deepParams, const MatrixConstView &data,
+           Matrix &Q, Matrix &H, size_t presampleStart, int &info, const VectorConstView &nruns, 
+           size_t fblock, size_t nBlocks, Proposal pdd, EstimatedParametersDescription &epd,
            const std::string &resultsFileStem, size_t console_mode, size_t load_mh_file)
 {
   enum {iMin, iMax};
@@ -382,7 +382,7 @@ sampleMHMC(LogPosteriorDensity &lpd, RandomWalkMetropolisHastings &rwmh,
           try
             {
               jsux = rwmh.compute(mhLogPostDens, mhParamDraws, steadyState, startParams, deepParams, data, Q, H,
-                              presampleStart, info, irun, currInitSizeArray, Jscale, lpd, drawDistribution, epd);
+                              presampleStart, info, irun, currInitSizeArray, lpd, pdd, epd);
               irun = currInitSizeArray;
               sux += jsux*currInitSizeArray;
               j += currInitSizeArray; //j=j+1;
@@ -715,18 +715,13 @@ logMCMCposterior(const VectorConstView &estParams, const MatrixConstView &data,
   GaussianPrior drawGaussDist01(0.0, 1.0, -INFINITY, INFINITY, 0.0, 1.0);
   // get Jscale = diag(bayestopt_.jscale);
   const mxArray *bayestopt_ = mexGetVariablePtr("global", "bayestopt_");
-  Matrix Jscale(n_estParams);
-  Matrix Dscale(n_estParams);
-  //Vector vJscale(n_estParams);
-  Jscale.setAll(0.0);
-  VectorConstView vJscale(mxGetPr(mxGetField(bayestopt_, 0, "jscale")), n_estParams, 1);
-  for (size_t i = 0; i < n_estParams; i++)
-    Jscale(i, i) = vJscale(i);
-  blas::gemm("N", "N", 1.0, D, Jscale, 0.0, Dscale);
+  const Matrix Jscale(n_estParams);
+  const VectorConstView vJscale(mxGetPr(mxGetField(bayestopt_, 0, "jscale")), n_estParams, 1);
+  Proposal pdd(vJscale,D);
 
   //sample MHMCMC draws and get get last line run in the last MH block sub-array
   int lastMHblockArrayLine = sampleMHMC(lpd, rwmh, steadyState, estParams2, deepParams, data, Q, H, presample, info,
-                                           nMHruns, fblock, nBlocks, Dscale, drawGaussDist01, epd, resultsFileStem, console_mode, load_mh_file);
+                                           nMHruns, fblock, nBlocks, pdd, epd, resultsFileStem, console_mode, load_mh_file);
 
   // Cleanups
   for (std::vector<EstimatedParameter>::iterator it = estParamsInfo.begin();

mex/sources/estimation/tests/random_walk_metropolis_hastings_core.m

@@ -136,7 +136,7 @@ oldoptions_console_mode=options_.console_mode;
     record.Seeds(b).Unifor = rand('state');
   end
   % calculate draws and get last line run in the last MH block sub-array
-  irun = logMHMCMCposterior( xparam1, varargin{2}, mexext, fblck, nblck, nruns, d)
+  irun = logMHMCMCposterior( xparam1, varargin{2}, mexext, fblck, nblck, nruns, vv)
   if irun<0
     error('Error in logMHMCMCposterior');
   end