stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

C++ Estimation DLL: Adding 1st cut draft RandomWalkMetropolisHastings.cc and related untested files for review

time-shift
George Perendia 2010-09-14 12:21:45 +01:00
parent f6aace084b
commit 99a884e09c
4 changed files with 214 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

43
mex/sources/estimation/RandSampler.cc Normal file
View File

@ -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 <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));
}

52
mex/sources/estimation/RandSampler.hh Normal file
View File

@ -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 <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_)

65
mex/sources/estimation/RandomWalkMetropolisHastings.cc Normal file
View File

@ -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 <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)
{
}

54
mex/sources/estimation/RandomWalkMetropolisHastings.hh Normal file
View File

@ -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 <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_)