From 99a884e09c26e5e13de0661018e03dc5254aa59e Mon Sep 17 00:00:00 2001 From: George Perendia Date: Tue, 14 Sep 2010 12:21:45 +0100 Subject: [PATCH] C++ Estimation DLL: Adding 1st cut draft RandomWalkMetropolisHastings.cc and related untested files for review --- mex/sources/estimation/RandSampler.cc | 43 ++++++++++++ mex/sources/estimation/RandSampler.hh | 52 +++++++++++++++ .../RandomWalkMetropolisHastings.cc | 65 +++++++++++++++++++ .../RandomWalkMetropolisHastings.hh | 54 +++++++++++++++ 4 files changed, 214 insertions(+) create mode 100644 mex/sources/estimation/RandSampler.cc create mode 100644 mex/sources/estimation/RandSampler.hh create mode 100644 mex/sources/estimation/RandomWalkMetropolisHastings.cc create mode 100644 mex/sources/estimation/RandomWalkMetropolisHastings.hh diff --git a/mex/sources/estimation/RandSampler.cc b/mex/sources/estimation/RandSampler.cc new file mode 100644 index 000000000..ddf904fe5 --- /dev/null +++ b/mex/sources/estimation/RandSampler.cc @@ -0,0 +1,43 @@ +/* + * 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 . + */ + +/////////////////////////////////////////////////////////// +// 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)); +} + diff --git a/mex/sources/estimation/RandSampler.hh b/mex/sources/estimation/RandSampler.hh new file mode 100644 index 000000000..16a08e427 --- /dev/null +++ b/mex/sources/estimation/RandSampler.hh @@ -0,0 +1,52 @@ +/* + * 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 . + */ + +/////////////////////////////////////////////////////////// +// 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_) diff --git a/mex/sources/estimation/RandomWalkMetropolisHastings.cc b/mex/sources/estimation/RandomWalkMetropolisHastings.cc new file mode 100644 index 000000000..b92bb268b --- /dev/null +++ b/mex/sources/estimation/RandomWalkMetropolisHastings.cc @@ -0,0 +1,65 @@ +/* + * 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 . + */ + +#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) +{ + +} + diff --git a/mex/sources/estimation/RandomWalkMetropolisHastings.hh b/mex/sources/estimation/RandomWalkMetropolisHastings.hh new file mode 100644 index 000000000..7507f4ef8 --- /dev/null +++ b/mex/sources/estimation/RandomWalkMetropolisHastings.hh @@ -0,0 +1,54 @@ +/* + * 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 . + */ + +/////////////////////////////////////////////////////////// +// 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_)