Merge remote branch 'george/master'

mex/sources/estimation/RandSampler.cc

+/*
+ * 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.cpp
+//  Implementation of the Class RandSampler
+//  Created on:      09-Sep-2010 23:25:30
+///////////////////////////////////////////////////////////
+
+#include "RandSampler.hh"
+
+/**
+ * draw = Mean + randn(1,n) * Sigma_upper_chol;
+ */
+void
+RandSampler::randMultiVar(Prior &distribution, Vector &draw, const Vector &mean, const Matrix &Scale, const size_t n)
+{
+  assert(n == draw.getSize());
+  assert(n == mean.getSize() || 1 == mean.getSize());
+  assert(n == Scale.getRows() || 1 == Scale.getRows());
+  assert(n == Scale.getCols() || 1 == Scale.getCols());
+  for (size_t i = 0; i < n; ++i)
+    draw(i) = mean(1 == mean.getSize() ? 0 : i)
+              + distribution.drand()
+              *Scale((1 == Scale.getRows() ? 0 : i), (1 == Scale.getCols() ? 0 : i));
+}
+

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:
+  RandSampler(){
+  };
+  virtual ~RandSampler(){};
+  virtual double compute(Vector &mhLogPostDens, MatrixView &mhParams, Matrix &steadyState,
+                         Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H,
+                         size_t presampleStart, int &info, const size_t nMHruns = 0, const Matrix &Jscale,
+                         const Matrix &D, const LogPosteriorDensity &logPosteriorDensity, const Prior &drawDistribution) = 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);
+
+  virtual void saveDraws() = 0;
+
+};
+
+#endif // !defined(5D7E5E52_2A4F_4f98_9A5C_A7FD8C278E0A__INCLUDED_)

mex/sources/estimation/RandomWalkMetropolisHastings.cc

+/*
+ * 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/>.
+ */
+
+#include "RandomWalkMetropolisHastings.hh"
+
+double
+RandomWalkMetropolisHastings::compute(Vector &mhLogPostDens, MatrixView &mhParams, Matrix &steadyState,
+                                      Vector &estParams2, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H,
+                                      const size_t presampleStart, int &info, const size_t nMHruns, const Matrix &Jscale,
+                                      const Matrix &D, LogPosteriorDensity &lpd, Prior &drawDistribution)
+{
+  double logpost, newLogpost;
+  size_t accepted = 0;
+  parDraw = estParams2;
+  blas::gemm("N", "N", 1.0, D, Jscale, 1.0, Dscale);
+  for (size_t run = 0; run < nMHruns; ++run)
+    {
+      randMultiVar(drawDistribution, newParDraw, parDraw, Dscale, parDraw.getSize());
+      try
+        {
+          newLogpost = lpd.compute(steadyState, newParDraw, deepParams, data, Q, H, presampleStart, info);
+        }
+      catch(...)
+        {
+          newLogpost = -INFINITY;
+        }
+      if ((newLogpost > -INFINITY) && log(uniform.drand()) < newLogpost-logpost)
+        {
+          mat::get_col(mhParams, run) = newParDraw;
+          parDraw = newParDraw;
+          mhLogPostDens(run) = newLogpost;
+          logpost = newLogpost;
+          accepted++;
+        }
+      else
+        {
+          mat::get_col(mhParams, run) = parDraw;
+          mhLogPostDens(run) = logpost;
+        }
+    }
+  return accepted/nMHruns;
+}
+
+void
+RandomWalkMetropolisHastings::saveDraws(const std::string &modName, const std::string &suffix, const MatrixView &Draws, const size_t block)
+{
+
+}
+

mex/sources/estimation/RandomWalkMetropolisHastings.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/>.
+ */
+
+///////////////////////////////////////////////////////////
+//  RandomWalkMetropolisHastings.hh
+//  Implementation of the Class RandomWalkMetropolisHastings
+//  Created on:      07-Sep-2010 15:21:40
+///////////////////////////////////////////////////////////
+
+#if !defined(A6BBC5E0_598E_4863_B7FF_E87320056B80__INCLUDED_)
+#define A6BBC5E0_598E_4863_B7FF_E87320056B80__INCLUDED_
+
+#include "RandSampler.hh"
+
+class RandomWalkMetropolisHastings : public RandSampler {
+
+private:
+  UniformPrior uniform;
+  Matrix Dscale;
+  Vector parDraw, newParDraw;
+
+public:
+  RandomWalkMetropolisHastings(size_t size) :
+    uniform(0.0, 0.0, 0.0, 1.0, 0.0, 1.0),
+    Dscale(size), parDraw(size), newParDraw(size)
+  {
+  };
+  virtual ~RandomWalkMetropolisHastings(){};
+  double compute(Vector &mhLogPostDens, MatrixView &mhParams, Matrix &steadyState,
+                 Vector &estParams2, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H,
+                 const size_t presampleStart, int &info, const size_t nMHruns, const Matrix &Jscale,
+                 const Matrix &D, LogPosteriorDensity &logPosteriorDensity, Prior &drawDistribution);
+//	void randMultiVar(const Prior &distribution, Vector &draw, const Matrix &Scale, const Vector &mean, const size_t n = 0);
+  void saveDraws(const std::string &modName, const std::string &suffix, const MatrixView &Draws, const size_t block);
+
+};
+
+#endif // !defined(A6BBC5E0_598E_4863_B7FF_E87320056B80__INCLUDED_)