From f89edf476e0a5c69fca0cad0236dda343717c7ca Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Wed, 6 Mar 2013 17:55:15 +0100 Subject: [PATCH] Start refactoring of penalty handling Handle all events that lead to a penalty through an exception. Next steps: - make all exceptions derive from a common class containing info code and penalty value - catch these exceptions at the top-level, and construct the penalty-adjusted likelihood --- .../estimation/InitializeKalmanFilter.cc | 31 ++--------- .../estimation/InitializeKalmanFilter.hh | 24 ++++----- mex/sources/estimation/KalmanFilter.cc | 27 +++++----- mex/sources/estimation/KalmanFilter.hh | 30 +++-------- mex/sources/estimation/LogLikelihoodMain.cc | 4 +- mex/sources/estimation/LogLikelihoodMain.hh | 6 +-- .../estimation/LogLikelihoodSubSample.cc | 6 +-- .../estimation/LogLikelihoodSubSample.hh | 54 +++++++------------ mex/sources/estimation/LogPosteriorDensity.cc | 4 +- mex/sources/estimation/LogPosteriorDensity.hh | 6 +-- .../RandomWalkMetropolisHastings.hh | 8 +-- mex/sources/estimation/logMHMCMCposterior.cc | 8 +-- mex/sources/estimation/logposterior.cc | 13 +++-- 13 files changed, 84 insertions(+), 137 deletions(-) diff --git a/mex/sources/estimation/InitializeKalmanFilter.cc b/mex/sources/estimation/InitializeKalmanFilter.cc index 0f989d70f..53ea5df6c 100644 --- a/mex/sources/estimation/InitializeKalmanFilter.cc +++ b/mex/sources/estimation/InitializeKalmanFilter.cc @@ -37,8 +37,7 @@ InitializeKalmanFilter::InitializeKalmanFilter(const std::string &basename, size const std::vector &varobs_arg, double qz_criterium_arg, double lyapunov_tol_arg, - bool noconstant_arg, - int &info) : + bool noconstant_arg) : lyapunov_tol(lyapunov_tol_arg), zeta_varobs_back_mixed(zeta_varobs_back_mixed_arg), detrendData(varobs_arg, noconstant_arg), @@ -60,7 +59,7 @@ InitializeKalmanFilter::InitializeKalmanFilter(const std::string &basename, size } void -InitializeKalmanFilter::setT(Matrix &T, int &info) +InitializeKalmanFilter::setT(Matrix &T) { // Initialize the empty columns of T to zero T.setAll(0.0); @@ -68,31 +67,11 @@ InitializeKalmanFilter::setT(Matrix &T, int &info) } void -InitializeKalmanFilter::setPstar(Matrix &Pstar, Matrix &Pinf, const Matrix &T, const Matrix &RQRt, int &info) +InitializeKalmanFilter::setPstar(Matrix &Pstar, Matrix &Pinf, const Matrix &T, const Matrix &RQRt) throw (DiscLyapFast::DLPException) { - try - { - // disclyap_fast(T, RQR, Pstar, lyapunov_tol, 0 or 1 to check chol) - discLyapFast.solve_lyap(T, RQRt, Pstar, lyapunov_tol, 0); + discLyapFast.solve_lyap(T, RQRt, Pstar, lyapunov_tol, 0); - Pinf.setAll(0.0); - } - catch (const DiscLyapFast::DLPException &e) - { - if (e.info > 0) // The matrix is not positive definite in NormCholesky calculator - { - puts(e.message.c_str()); - info = -1; //likelihood = penalty; - return; - } - else if (e.info < 0) - { - printf("Caugth unhandled TS exception with Pstar matrix: "); - puts(e.message.c_str()); - info = -1; //likelihood = penalty; - throw; - } - } + Pinf.setAll(0.0); } diff --git a/mex/sources/estimation/InitializeKalmanFilter.hh b/mex/sources/estimation/InitializeKalmanFilter.hh index de266d2d1..83e25a19f 100644 --- a/mex/sources/estimation/InitializeKalmanFilter.hh +++ b/mex/sources/estimation/InitializeKalmanFilter.hh @@ -48,31 +48,31 @@ public: const std::vector &zeta_varobs_back_mixed_arg, const std::vector &varobs_arg, double qz_criterium_arg, double lyapunov_tol_arg, - bool noconstant_arg, int &info); + bool noconstant_arg); virtual ~InitializeKalmanFilter(); // initialise parameter dependent KF matrices only but not Ps template void initialize(Vec1 &steadyState, const Vec2 &deepParams, Mat1 &R, const Mat2 &Q, Matrix &RQRt, Matrix &T, - double &penalty, const MatrixConstView &dataView, - MatrixView &detrendedDataView, int &info) + const MatrixConstView &dataView, + MatrixView &detrendedDataView) { modelSolution.compute(steadyState, deepParams, g_x, g_u); detrendData.detrend(steadyState, dataView, detrendedDataView); - setT(T, info); - setRQR(R, Q, RQRt, info); + setT(T); + setRQR(R, Q, RQRt); } // initialise all KF matrices template void initialize(Vec1 &steadyState, const Vec2 &deepParams, Mat1 &R, const Mat2 &Q, Matrix &RQRt, Matrix &T, Matrix &Pstar, Matrix &Pinf, - double &penalty, const MatrixConstView &dataView, - MatrixView &detrendedDataView, int &info) + const MatrixConstView &dataView, + MatrixView &detrendedDataView) { - initialize(steadyState, deepParams, R, Q, RQRt, T, penalty, dataView, detrendedDataView, info); - setPstar(Pstar, Pinf, T, RQRt, info); + initialize(steadyState, deepParams, R, Q, RQRt, T, dataView, detrendedDataView); + setPstar(Pstar, Pinf, T, RQRt); } private: @@ -87,10 +87,10 @@ private: Matrix g_x; Matrix g_u; Matrix Rt, RQ; - void setT(Matrix &T, int &info); + void setT(Matrix &T); template - void setRQR(Mat1 &R, const Mat2 &Q, Matrix &RQRt, int &info) + void setRQR(Mat1 &R, const Mat2 &Q, Matrix &RQRt) { mat::assignByVectors(R, mat::nullVec, mat::nullVec, g_u, zeta_varobs_back_mixed, mat::nullVec); @@ -98,7 +98,7 @@ private: blas::gemm("N", "N", 1.0, R, Q, 0.0, RQ); // R*Q blas::gemm("N", "T", 1.0, RQ, R, 0.0, RQRt); // R*Q*R' } - void setPstar(Matrix &Pstar, Matrix &Pinf, const Matrix &T, const Matrix &RQRt, int &info); + void setPstar(Matrix &Pstar, Matrix &Pinf, const Matrix &T, const Matrix &RQRt) throw (DiscLyapFast::DLPException); }; diff --git a/mex/sources/estimation/KalmanFilter.cc b/mex/sources/estimation/KalmanFilter.cc index b624da675..32f121899 100644 --- a/mex/sources/estimation/KalmanFilter.cc +++ b/mex/sources/estimation/KalmanFilter.cc @@ -1,5 +1,5 @@ /* - * Copyright (C) 2009-2012 Dynare Team + * Copyright (C) 2009-2013 Dynare Team * * This file is part of Dynare. * @@ -36,7 +36,7 @@ KalmanFilter::KalmanFilter(const std::string &basename, size_t n_endo, size_t n_ const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, double qz_criterium_arg, const std::vector &varobs_arg, double riccati_tol_arg, double lyapunov_tol_arg, - bool noconstant_arg, int &info) : + bool noconstant_arg) : zeta_varobs_back_mixed(compute_zeta_varobs_back_mixed(zeta_back_arg, zeta_mixed_arg, varobs_arg)), Z(varobs_arg.size(), zeta_varobs_back_mixed.size()), Zt(Z.getCols(), Z.getRows()), T(zeta_varobs_back_mixed.size()), R(zeta_varobs_back_mixed.size(), n_exo), Pstar(zeta_varobs_back_mixed.size(), zeta_varobs_back_mixed.size()), Pinf(zeta_varobs_back_mixed.size(), zeta_varobs_back_mixed.size()), @@ -45,7 +45,7 @@ KalmanFilter::KalmanFilter(const std::string &basename, size_t n_endo, size_t n_ oldKFinv(zeta_varobs_back_mixed.size(), varobs_arg.size()), a_init(zeta_varobs_back_mixed.size()), a_new(zeta_varobs_back_mixed.size()), vt(varobs_arg.size()), vtFinv(varobs_arg.size()), riccati_tol(riccati_tol_arg), initKalmanFilter(basename, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, - zeta_static_arg, zeta_varobs_back_mixed, varobs_arg, qz_criterium_arg, lyapunov_tol_arg, noconstant_arg, info), + zeta_static_arg, zeta_varobs_back_mixed, varobs_arg, qz_criterium_arg, lyapunov_tol_arg, noconstant_arg), FUTP(varobs_arg.size()*(varobs_arg.size()+1)/2) { Z.setAll(0.0); @@ -79,12 +79,14 @@ KalmanFilter::compute_zeta_varobs_back_mixed(const std::vector &zeta_bac * Multi-variate standard Kalman Filter */ double -KalmanFilter::filter(const MatrixView &detrendedDataView, const Matrix &H, VectorView &vll, size_t start, int &info) +KalmanFilter::filter(const MatrixView &detrendedDataView, const Matrix &H, VectorView &vll, size_t start) { double loglik = 0.0, ll, logFdet = 0.0, Fdet, dvtFinvVt; size_t p = Finv.getRows(); bool nonstationary = true; a_init.setAll(0.0); + int info; + for (size_t t = 0; t < detrendedDataView.getCols(); ++t) { if (nonstationary) @@ -106,12 +108,8 @@ KalmanFilter::filter(const MatrixView &detrendedDataView, const Matrix &H, Vect FUTP(i + (j-1)*j/2 -1) = F(i-1, j-1); info = lapack::choleskySolver(FUTP, Finv, "U"); // F now contains its Chol decomposition! - if (info < 0) - { - std::cout << "Info:" << info << std::endl; - std::cout << "t:" << t << std::endl; - return 0; - } + assert(info >= 0); + if (info > 0) { //enforce Pstar symmetry with P=(P+P')/2=0.5P+0.5P' @@ -136,11 +134,10 @@ KalmanFilter::filter(const MatrixView &detrendedDataView, const Matrix &H, Vect for (size_t j = i; j <= p; ++j) FUTP(i + (j-1)*j/2 -1) = F(i-1, j-1); - info = lapack::choleskySolver(FUTP, Finv, "U"); // F now contains its Chol decomposition! - if (info != 0) - { - return 0; - } + info = lapack::choleskySolver(FUTP, Finv, "U"); // F now contains + // its Chol + // decomposition! + assert(info == 0); } // KFinv gain matrix blas::symm("R", "U", 1.0, Finv, K, 0.0, KFinv); diff --git a/mex/sources/estimation/KalmanFilter.hh b/mex/sources/estimation/KalmanFilter.hh index 4184efe78..b9203d610 100644 --- a/mex/sources/estimation/KalmanFilter.hh +++ b/mex/sources/estimation/KalmanFilter.hh @@ -53,38 +53,24 @@ public: const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, double qz_criterium_arg, const std::vector &varobs_arg, double riccati_tol_arg, double lyapunov_tol_arg, - bool noconstant_arg, int &info); + bool noconstant_arg); template double compute(const MatrixConstView &dataView, Vec1 &steadyState, const Mat1 &Q, const Matrix &H, const Vec2 &deepParams, - VectorView &vll, MatrixView &detrendedDataView, size_t start, size_t period, - double &penalty, int &info) + VectorView &vll, MatrixView &detrendedDataView, size_t start, size_t period) { double lik = INFINITY; - try - { + if (period == 0) // initialise all KF matrices initKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T, Pstar, Pinf, - penalty, dataView, detrendedDataView, info); + dataView, detrendedDataView); else // initialise parameter dependent KF matrices only but not Ps initKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T, - penalty, dataView, detrendedDataView, info); + dataView, detrendedDataView); - lik = filter(detrendedDataView, H, vll, start, info); - } - catch (const DecisionRules::BlanchardKahnException &bke) - { - info = 22; - return penalty; - } - - if (info != 0) - return penalty; - else - return lik; - - }; + return filter(detrendedDataView, H, vll, start); + } private: const std::vector zeta_varobs_back_mixed; @@ -106,7 +92,7 @@ private: Vector FUTP; // F upper triangle packed as vector FUTP(i + (j-1)*j/2) = F(i,j) for 1<=i<=j; // Method - double filter(const MatrixView &detrendedDataView, const Matrix &H, VectorView &vll, size_t start, int &info); + double filter(const MatrixView &detrendedDataView, const Matrix &H, VectorView &vll, size_t start); }; diff --git a/mex/sources/estimation/LogLikelihoodMain.cc b/mex/sources/estimation/LogLikelihoodMain.cc index deaff883b..8109ad99a 100644 --- a/mex/sources/estimation/LogLikelihoodMain.cc +++ b/mex/sources/estimation/LogLikelihoodMain.cc @@ -29,11 +29,11 @@ LogLikelihoodMain::LogLikelihoodMain(const std::string &basename, EstimatedParam const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, const double qz_criterium, const std::vector &varobs, double riccati_tol, double lyapunov_tol, - bool noconstant_arg, int &info_arg) + bool noconstant_arg) : estSubsamples(estiParDesc.estSubsamples), logLikelihoodSubSample(basename, estiParDesc, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, qz_criterium, - varobs, riccati_tol, lyapunov_tol, noconstant_arg, info_arg), + varobs, riccati_tol, lyapunov_tol, noconstant_arg), vll(estiParDesc.getNumberOfPeriods()), // time dimension size of data detrendedData(varobs.size(), estiParDesc.getNumberOfPeriods()) { diff --git a/mex/sources/estimation/LogLikelihoodMain.hh b/mex/sources/estimation/LogLikelihoodMain.hh index a2f062888..87296e147 100644 --- a/mex/sources/estimation/LogLikelihoodMain.hh +++ b/mex/sources/estimation/LogLikelihoodMain.hh @@ -41,7 +41,7 @@ public: const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, const double qz_criterium_arg, const std::vector &varobs_arg, double riccati_tol_arg, double lyapunov_tol_arg, - bool noconstant_arg, int &info); + bool noconstant_arg); /** * Compute method Inputs: @@ -54,7 +54,7 @@ public: template double compute(VEC1 &steadyState, VEC2 &estParams, VectorView &deepParams, const MatrixConstView &data, - MatrixView &Q, Matrix &H, size_t start, int &info) + MatrixView &Q, Matrix &H, size_t start) { double logLikelihood = 0; for (size_t i = 0; i < estSubsamples.size(); ++i) @@ -68,7 +68,7 @@ public: VectorView vllView(vll, estSubsamples[i].startPeriod, estSubsamples[i].endPeriod-estSubsamples[i].startPeriod+1); logLikelihood += logLikelihoodSubSample.compute(vSteadyState, dataView, estParams, deepParams, - Q, H, vllView, detrendedDataView, info, start, i); + Q, H, vllView, detrendedDataView, start, i); } return logLikelihood; }; diff --git a/mex/sources/estimation/LogLikelihoodSubSample.cc b/mex/sources/estimation/LogLikelihoodSubSample.cc index 0a391b1bc..9a3aabf80 100644 --- a/mex/sources/estimation/LogLikelihoodSubSample.cc +++ b/mex/sources/estimation/LogLikelihoodSubSample.cc @@ -32,10 +32,10 @@ LogLikelihoodSubSample::~LogLikelihoodSubSample() LogLikelihoodSubSample::LogLikelihoodSubSample(const std::string &basename, EstimatedParametersDescription &INestiParDesc, size_t n_endo, size_t n_exo, const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, const double qz_criterium, - const std::vector &varobs, double riccati_tol, double lyapunov_tol, bool noconstant_arg, int &INinfo) : - startPenalty(-1e8), estiParDesc(INestiParDesc), + const std::vector &varobs, double riccati_tol, double lyapunov_tol, bool noconstant_arg) : + estiParDesc(INestiParDesc), kalmanFilter(basename, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, qz_criterium, - varobs, riccati_tol, lyapunov_tol, noconstant_arg, INinfo), eigQ(n_exo), eigH(varobs.size()), info(INinfo) + varobs, riccati_tol, lyapunov_tol, noconstant_arg), eigQ(n_exo), eigH(varobs.size()) { }; diff --git a/mex/sources/estimation/LogLikelihoodSubSample.hh b/mex/sources/estimation/LogLikelihoodSubSample.hh index 8e8502be8..4e0fa509f 100644 --- a/mex/sources/estimation/LogLikelihoodSubSample.hh +++ b/mex/sources/estimation/LogLikelihoodSubSample.hh @@ -39,35 +39,33 @@ public: LogLikelihoodSubSample(const std::string &basename, EstimatedParametersDescription &estiParDesc, size_t n_endo, size_t n_exo, const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, const double qz_criterium, - const std::vector &varobs_arg, double riccati_tol_in, double lyapunov_tol, bool noconstant_arg, int &info); + const std::vector &varobs_arg, double riccati_tol_in, double lyapunov_tol, bool noconstant_arg); template double compute(VEC1 &steadyState, const MatrixConstView &dataView, VEC2 &estParams, VectorView &deepParams, - MatrixView &Q, Matrix &H, VectorView &vll, MatrixView &detrendedDataView, int &info, size_t start, size_t period) + MatrixView &Q, Matrix &H, VectorView &vll, MatrixView &detrendedDataView, size_t start, size_t period) { - penalty = startPenalty; - logLikelihood = startPenalty; - updateParams(estParams, deepParams, Q, H, period); - if (info == 0) - logLikelihood = kalmanFilter.compute(dataView, steadyState, Q, H, deepParams, vll, detrendedDataView, start, period, penalty, info); - // else - // logLikelihood+=penalty; - return logLikelihood; - - }; + return kalmanFilter.compute(dataView, steadyState, Q, H, deepParams, vll, detrendedDataView, start, period); + } virtual ~LogLikelihoodSubSample(); + class UpdateParamsException + { + public: + double penalty; + UpdateParamsException(double penalty_arg) : penalty(penalty_arg) + { + } + }; + private: - double startPenalty, penalty; - double logLikelihood; EstimatedParametersDescription &estiParDesc; KalmanFilter kalmanFilter; VDVEigDecomposition eigQ; VDVEigDecomposition eigH; - int &info; // methods template @@ -78,7 +76,6 @@ private: int test; bool found; std::vector::const_iterator it; - info = 0; for (i = 0; i < estParams.getSize(); ++i) { @@ -108,13 +105,9 @@ private: Q(k2, k1) = Q(k1, k2); // [CholQ,testQ] = chol(Q); test = lapack::choleskyDecomp(Q, "L"); + assert(test >= 0); + if (test > 0) - { - mexPrintf("Caught unhandled exception with cholesky of Q matrix: "); - logLikelihood = penalty; - info = 1; - } - else if (test < 0) { // The variance-covariance matrix of the structural innovations is not definite positive. // We have to compute the eigenvalues of this matrix in order to build the penalty. @@ -129,8 +122,7 @@ private: delta -= evQ(i); } - logLikelihood = penalty+delta; - info = 43; + throw UpdateParamsException(delta); } // if break; @@ -144,13 +136,9 @@ private: //[CholH,testH] = chol(H); test = lapack::choleskyDecomp(H, "L"); + assert(test >= 0); + if (test > 0) - { - mexPrintf("Caught unhandled exception with cholesky of Q matrix: "); - logLikelihood = penalty; - info = 1; - } - else if (test < 0) { // The variance-covariance matrix of the measurement errors is not definite positive. // We have to compute the eigenvalues of this matrix in order to build the penalty. @@ -165,8 +153,7 @@ private: if (evH(i) < 0) delta -= evH(i); } - logLikelihood = penalty+delta; - info = 44; + throw UpdateParamsException(delta); } // end if break; @@ -176,8 +163,7 @@ private: deepParams(k) = estParams(i); break; default: - logLikelihood = penalty; - info = 1; + assert(false); } // end switch } // end found } //end for diff --git a/mex/sources/estimation/LogPosteriorDensity.cc b/mex/sources/estimation/LogPosteriorDensity.cc index 114be39ea..a9da75fc1 100644 --- a/mex/sources/estimation/LogPosteriorDensity.cc +++ b/mex/sources/estimation/LogPosteriorDensity.cc @@ -33,10 +33,10 @@ LogPosteriorDensity::LogPosteriorDensity(const std::string &modName, EstimatedPa const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, const double qz_criterium_arg, const std::vector &varobs_arg, double riccati_tol_arg, double lyapunov_tol_arg, - bool noconstant_arg, int &info_arg) : + bool noconstant_arg) : logPriorDensity(estParamsDesc), logLikelihoodMain(modName, estParamsDesc, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, - zeta_static_arg, qz_criterium_arg, varobs_arg, riccati_tol_arg, lyapunov_tol_arg, noconstant_arg, info_arg) + zeta_static_arg, qz_criterium_arg, varobs_arg, riccati_tol_arg, lyapunov_tol_arg, noconstant_arg) { } diff --git a/mex/sources/estimation/LogPosteriorDensity.hh b/mex/sources/estimation/LogPosteriorDensity.hh index e11795c9e..820f2a575 100644 --- a/mex/sources/estimation/LogPosteriorDensity.hh +++ b/mex/sources/estimation/LogPosteriorDensity.hh @@ -48,13 +48,13 @@ public: const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, const double qz_criterium_arg, const std::vector &varobs_arg, double riccati_tol_arg, double lyapunov_tol_arg, - bool noconstant_arg, int &info_arg); + bool noconstant_arg); template double - compute(VEC1 &steadyState, VEC2 &estParams, VectorView &deepParams, const MatrixConstView &data, MatrixView &Q, Matrix &H, size_t presampleStart, int &info) + compute(VEC1 &steadyState, VEC2 &estParams, VectorView &deepParams, const MatrixConstView &data, MatrixView &Q, Matrix &H, size_t presampleStart) { - return -logLikelihoodMain.compute(steadyState, estParams, deepParams, data, Q, H, presampleStart, info) + return -logLikelihoodMain.compute(steadyState, estParams, deepParams, data, Q, H, presampleStart) -logPriorDensity.compute(estParams); } diff --git a/mex/sources/estimation/RandomWalkMetropolisHastings.hh b/mex/sources/estimation/RandomWalkMetropolisHastings.hh index 5b5c02e9b..f281939ba 100644 --- a/mex/sources/estimation/RandomWalkMetropolisHastings.hh +++ b/mex/sources/estimation/RandomWalkMetropolisHastings.hh @@ -1,5 +1,5 @@ /* - * Copyright (C) 2010-2012 Dynare Team + * Copyright (C) 2010-2013 Dynare Team * * This file is part of Dynare. * @@ -47,7 +47,7 @@ public: template double compute(VectorView &mhLogPostDens, MatrixView &mhParams, VEC1 &steadyState, Vector &estParams, VectorView &deepParams, const MatrixConstView &data, MatrixView &Q, Matrix &H, - const size_t presampleStart, int &info, const size_t startDraw, size_t nMHruns, + const size_t presampleStart, const size_t startDraw, size_t nMHruns, LogPosteriorDensity &lpd, Proposal &pDD, EstimatedParametersDescription &epd) { //streambuf *likbuf, *drawbuf *backup; @@ -60,7 +60,7 @@ public: size_t count, accepted = 0; parDraw = estParams; - logpost = -lpd.compute(steadyState, estParams, deepParams, data, Q, H, presampleStart, info); + logpost = -lpd.compute(steadyState, estParams, deepParams, data, Q, H, presampleStart); for (size_t run = startDraw - 1; run < nMHruns; ++run) { @@ -79,7 +79,7 @@ public: { try { - newLogpost = -lpd.compute(steadyState, newParDraw, deepParams, data, Q, H, presampleStart, info); + newLogpost = -lpd.compute(steadyState, newParDraw, deepParams, data, Q, H, presampleStart); } catch (const std::exception &e) { diff --git a/mex/sources/estimation/logMHMCMCposterior.cc b/mex/sources/estimation/logMHMCMCposterior.cc index 1d4225d22..68a519b61 100644 --- a/mex/sources/estimation/logMHMCMCposterior.cc +++ b/mex/sources/estimation/logMHMCMCposterior.cc @@ -114,7 +114,7 @@ fillEstParamsInfo(const mxArray *estim_params_info, EstimatedParameter::pType ty int sampleMHMC(LogPosteriorDensity &lpd, RandomWalkMetropolisHastings &rwmh, VectorView &steadyState, VectorConstView &estParams, VectorView &deepParams, const MatrixConstView &data, - MatrixView &Q, Matrix &H, size_t presampleStart, int &info, const VectorConstView &nruns, + MatrixView &Q, Matrix &H, size_t presampleStart, 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) { @@ -397,7 +397,7 @@ sampleMHMC(LogPosteriorDensity &lpd, RandomWalkMetropolisHastings &rwmh, try { jsux = rwmh.compute(mhLogPostDens, mhParamDraws, steadyState, startParams, deepParams, data, Q, H, - presampleStart, info, irun, currInitSizeArray, lpd, pdd, epd); + presampleStart, irun, currInitSizeArray, lpd, pdd, epd); irun = currInitSizeArray; sux += jsux*currInitSizeArray; j += currInitSizeArray; //j=j+1; @@ -715,7 +715,7 @@ logMCMCposterior(VectorConstView &estParams, const MatrixConstView &data, // Allocate LogPosteriorDensity object int info; LogPosteriorDensity lpd(basename, epd, n_endo, n_exo, zeta_fwrd, zeta_back, zeta_mixed, zeta_static, - qz_criterium, varobs, riccati_tol, lyapunov_tol, noconstant, info); + qz_criterium, varobs, riccati_tol, lyapunov_tol, noconstant); // Construct MHMCMC Sampler RandomWalkMetropolisHastings rwmh(estParams.getSize()); @@ -727,7 +727,7 @@ logMCMCposterior(VectorConstView &estParams, const MatrixConstView &data, 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, estParams, deepParams, data, Q, H, presample, info, + int lastMHblockArrayLine = sampleMHMC(lpd, rwmh, steadyState, estParams, deepParams, data, Q, H, presample, nMHruns, fblock, nBlocks, pdd, epd, resultsFileStem, console_mode, load_mh_file); // Cleanups diff --git a/mex/sources/estimation/logposterior.cc b/mex/sources/estimation/logposterior.cc index 7add29b96..4fad4e2f8 100644 --- a/mex/sources/estimation/logposterior.cc +++ b/mex/sources/estimation/logposterior.cc @@ -104,7 +104,7 @@ double logposterior(VEC1 &estParams, const MatrixConstView &data, const mxArray *options_, const mxArray *M_, const mxArray *estim_params_, const mxArray *bayestopt_, const mxArray *oo_, VEC2 &steadyState, double *trend_coeff, - int &info, VectorView &deepParams, Matrix &H, MatrixView &Q) + VectorView &deepParams, Matrix &H, MatrixView &Q) { double loglinear = *mxGetPr(mxGetField(options_, 0, "loglinear")); if (loglinear == 1) @@ -183,12 +183,12 @@ logposterior(VEC1 &estParams, const MatrixConstView &data, // Allocate LogPosteriorDensity object LogPosteriorDensity lpd(basename, epd, n_endo, n_exo, zeta_fwrd, zeta_back, zeta_mixed, zeta_static, - qz_criterium, varobs, riccati_tol, lyapunov_tol, noconstant, info); + qz_criterium, varobs, riccati_tol, lyapunov_tol, noconstant); // Construct arguments of compute() method // Compute the posterior - double logPD = lpd.compute(steadyState, estParams, deepParams, data, Q, H, presample, info); + double logPD = lpd.compute(steadyState, estParams, deepParams, data, Q, H, presample); // Cleanups for (std::vector::iterator it = estParamsInfo.begin(); @@ -268,11 +268,10 @@ mexFunction(int nlhs, mxArray *plhs[], // Compute and return the value try { - int info; *lik = logposterior(estParams, data, options_, M_, estim_params_, bayestopt_, oo_, - steadyState, trend_coeff, info, deepParams, H, Q); - *info_mx = info; - *exit_flag = info; + steadyState, trend_coeff, deepParams, H, Q); + *info_mx = 0; + *exit_flag = 0; } catch (LogposteriorMexErrMsgTxtException e) {