From 4f98b560adec179c5d51e008613f5e6e16ab0e23 Mon Sep 17 00:00:00 2001 From: George Perendia Date: Sun, 30 May 2010 09:33:49 +0100 Subject: [PATCH] Estimation C++DLL: First draft for LogPosteriorDensity, LogPriorDensity and associated changes: random number generation still missing --- mex/sources/estimation/EstimatedParameter.cc | 2 +- mex/sources/estimation/EstimatedParameter.hh | 4 +- mex/sources/estimation/LogLikelihoodMain.hh | 3 +- .../estimation/LogLikelihoodSubSample.hh | 2 +- mex/sources/estimation/LogPosteriorDensity.cc | 67 +++++ mex/sources/estimation/LogPosteriorDensity.hh | 58 ++++ mex/sources/estimation/LogPriorDensity.cc | 58 ++++ mex/sources/estimation/LogPriorDensity.hh | 34 +++ mex/sources/estimation/Prior.cc | 4 +- mex/sources/estimation/Prior.hh | 259 +++++++++++++++--- 10 files changed, 439 insertions(+), 52 deletions(-) create mode 100644 mex/sources/estimation/LogPosteriorDensity.cc create mode 100644 mex/sources/estimation/LogPosteriorDensity.hh create mode 100644 mex/sources/estimation/LogPriorDensity.cc create mode 100644 mex/sources/estimation/LogPriorDensity.hh diff --git a/mex/sources/estimation/EstimatedParameter.cc b/mex/sources/estimation/EstimatedParameter.cc index b66b576f6..a61db1cb3 100644 --- a/mex/sources/estimation/EstimatedParameter.cc +++ b/mex/sources/estimation/EstimatedParameter.cc @@ -27,7 +27,7 @@ EstimatedParameter::EstimatedParameter(const EstimatedParameter::pType type_arg, size_t ID1_arg, size_t ID2_arg, const std::vector &subSampleIDs_arg, - double lower_bound_arg, double upper_bound_arg, Prior prior_arg) : + double lower_bound_arg, double upper_bound_arg, Prior* prior_arg) : ptype(type_arg), ID1(ID1_arg), ID2(ID2_arg), lower_bound(lower_bound_arg), upper_bound(upper_bound_arg), prior(prior_arg), subSampleIDs(subSampleIDs_arg) diff --git a/mex/sources/estimation/EstimatedParameter.hh b/mex/sources/estimation/EstimatedParameter.hh index 3a83b6270..2778220fd 100644 --- a/mex/sources/estimation/EstimatedParameter.hh +++ b/mex/sources/estimation/EstimatedParameter.hh @@ -43,7 +43,7 @@ public: EstimatedParameter(const EstimatedParameter::pType type, size_t ID1, size_t ID2, const std::vector &subSampleIDs, - double lower_bound, double upper_bound, Prior prior + double lower_bound, double upper_bound, Prior* prior ); virtual ~EstimatedParameter(); @@ -52,7 +52,7 @@ public: const size_t ID2; const double lower_bound; const double upper_bound; - const Prior prior; + Prior* prior; const std::vector subSampleIDs; }; diff --git a/mex/sources/estimation/LogLikelihoodMain.hh b/mex/sources/estimation/LogLikelihoodMain.hh index 2f428b4e6..cf790df7d 100644 --- a/mex/sources/estimation/LogLikelihoodMain.hh +++ b/mex/sources/estimation/LogLikelihoodMain.hh @@ -50,7 +50,8 @@ public: * Q and H KF matrices of shock and measurement error varinaces and covariances * KF logLikelihood calculation start period. */ - double compute(Matrix &steadyState, const Vector &estParams, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H, size_t presampleStart, int &info); // for calls from estimation and to set Steady State + double compute(Matrix &steadyState, const Vector &estParams, Vector &deepParams, const MatrixConstView &data, + Matrix &Q, Matrix &H, size_t presampleStart, int &info); // for calls from estimation and to set Steady State Vector & getVll() diff --git a/mex/sources/estimation/LogLikelihoodSubSample.hh b/mex/sources/estimation/LogLikelihoodSubSample.hh index 29a28f67a..bf92bfe28 100644 --- a/mex/sources/estimation/LogLikelihoodSubSample.hh +++ b/mex/sources/estimation/LogLikelihoodSubSample.hh @@ -26,8 +26,8 @@ #if !defined(DF8B7AF5_8169_4587_9037_2CD2C82E2DDF__INCLUDED_) #define DF8B7AF5_8169_4587_9037_2CD2C82E2DDF__INCLUDED_ -#include "KalmanFilter.hh" #include "EstimatedParametersDescription.hh" +#include "KalmanFilter.hh" #include "VDVEigDecomposition.hh" class LogLikelihoodSubSample { diff --git a/mex/sources/estimation/LogPosteriorDensity.cc b/mex/sources/estimation/LogPosteriorDensity.cc new file mode 100644 index 000000000..b3d2a204c --- /dev/null +++ b/mex/sources/estimation/LogPosteriorDensity.cc @@ -0,0 +1,67 @@ +/* + * 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 . + */ + +/////////////////////////////////////////////////////////// +// LogPosteriorDensity.cpp +// Implementation of the Class LogPosteriorDensity +// Created on: 10-Feb-2010 20:54:18 +/////////////////////////////////////////////////////////// + +#include "LogPosteriorDensity.hh" + +LogPosteriorDensity::~LogPosteriorDensity() +{ +} + +LogPosteriorDensity::LogPosteriorDensity(const std::string &modName, EstimatedParametersDescription &estParamsDesc, 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_arg, const std::vector &varobs_arg, + double riccati_tol_arg, double lyapunov_tol_arg, int &info_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, info_arg) +{ + +} + +double +LogPosteriorDensity::compute(Matrix &steadyState, const Vector &estParams, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H, size_t presampleStart, int &info) +{ + return logLikelihoodMain.compute(steadyState, estParams, deepParams, data, Q, H, presampleStart, info) + +logPriorDensity.compute(estParams); +} + +/** + * vector of log likelihoods for each Kalman step + */ +Vector & +LogPosteriorDensity::getLikVector() +{ + return logLikelihoodMain.getVll(); +} + +/** + * log likelihood as summ of Vll for each Kalman step + */ +double +LogPosteriorDensity::getLogPosteriorDensity() +{ + return logLikelihoodMain.getLogLikelihood()+logPriorDensity.getLogPriorDensity(); +} + diff --git a/mex/sources/estimation/LogPosteriorDensity.hh b/mex/sources/estimation/LogPosteriorDensity.hh new file mode 100644 index 000000000..fd4437919 --- /dev/null +++ b/mex/sources/estimation/LogPosteriorDensity.hh @@ -0,0 +1,58 @@ +/* + * 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 . + */ + +/////////////////////////////////////////////////////////// +// LogPosteriorDensity.hh +// Implementation of the Class LogPosteriorDensity +// Created on: 10-Feb-2010 20:54:18 +/////////////////////////////////////////////////////////// + +#if !defined(LPD_052A31B5_53BF_4904_AD80_863B52827973__INCLUDED_) +#define LPD_052A31B5_53BF_4904_AD80_863B52827973__INCLUDED_ + +#include "EstimatedParametersDescription.hh" +#include "LogPriorDensity.hh" +#include "LogLikelihoodMain.hh" + +/** + * Class that calculates Log Posterior Density using kalman, based on Dynare + * DsgeLikelihood.m + */ +class LogPosteriorDensity { + +private: + LogPriorDensity logPriorDensity; + LogLikelihoodMain logLikelihoodMain; + double logPosteriorDensity; + +public: + virtual ~LogPosteriorDensity(); + + LogPosteriorDensity(const std::string &modName, EstimatedParametersDescription &estParamsDesc, 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_arg, const std::vector &varobs_arg, + double riccati_tol_arg, double lyapunov_tol_arg, int &info_arg); + double compute(Matrix &steadyState, const Vector &estParams, Vector &deepParams, + const MatrixConstView &data, Matrix &Q, Matrix &H, size_t presampleStart, int &info); + Vector&getLikVector(); + double getLogPosteriorDensity(); + +}; + +#endif // !defined(052A31B5_53BF_4904_AD80_863B52827973__INCLUDED_) diff --git a/mex/sources/estimation/LogPriorDensity.cc b/mex/sources/estimation/LogPriorDensity.cc new file mode 100644 index 000000000..fd4577b87 --- /dev/null +++ b/mex/sources/estimation/LogPriorDensity.cc @@ -0,0 +1,58 @@ +/* + * 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 . + */ + +/////////////////////////////////////////////////////////// +// LogPriorDensity.cpp +// Implementation of the Class LogPriorDensity +// Created on: 10-Feb-2010 20:56:08 +/////////////////////////////////////////////////////////// + +#include "LogPriorDensity.hh" + +LogPriorDensity::~LogPriorDensity() +{ +}; + +LogPriorDensity::LogPriorDensity(EstimatedParametersDescription &estParsDesc_arg) : + estParsDesc(estParsDesc_arg) +{ +}; + +double +LogPriorDensity::compute(const Vector &ep) +{ + assert(estParsDesc.estParams.size() == ep.getSize()); + + for (size_t i = 0; i < ep.getSize(); ++i) + { + logPriorDensity += ((*(estParsDesc.estParams[i]).prior)).pdf(ep(i)); + if (std::isinf(abs(logPriorDensity))) + return logPriorDensity; + } + return logPriorDensity; +}; + +/** + * Return random number for prior fromits distribution + */ +void +LogPriorDensity::computeNewParams(Vector &ep) +{ + +} diff --git a/mex/sources/estimation/LogPriorDensity.hh b/mex/sources/estimation/LogPriorDensity.hh new file mode 100644 index 000000000..9732dd9d1 --- /dev/null +++ b/mex/sources/estimation/LogPriorDensity.hh @@ -0,0 +1,34 @@ +/////////////////////////////////////////////////////////// +// LogPriorDensity.hh +// Implementation of the Class LogPriorDensity +// Created on: 10-Feb-2010 20:56:08 +/////////////////////////////////////////////////////////// + +#if !defined(LPD_011FD4CF_17CE_4805_882B_046AA07CF443__INCLUDED_) +#define LPD_011FD4CF_17CE_4805_882B_046AA07CF443__INCLUDED_ + +//#include +#include "Vector.hh" +#include "EstimatedParametersDescription.hh" + +class LogPriorDensity { + +public: + LogPriorDensity(EstimatedParametersDescription& estParsDesc); + virtual ~LogPriorDensity(); + + double compute(const Vector &estParams); + double + getLogPriorDensity() + { + return logPriorDensity; + }; + void computeNewParams(Vector &newParams); + +private: + double logPriorDensity; + const EstimatedParametersDescription &estParsDesc; + +}; + +#endif // !defined(011FD4CF_17CE_4805_882B_046AA07CF443__INCLUDED_) diff --git a/mex/sources/estimation/Prior.cc b/mex/sources/estimation/Prior.cc index 0ca4dc58b..437c2ef04 100644 --- a/mex/sources/estimation/Prior.cc +++ b/mex/sources/estimation/Prior.cc @@ -25,8 +25,8 @@ #include "Prior.hh" -Prior::Prior(Prior::pShape shape_arg, double mean_arg, double mode_arg, double standard_arg, double lower_bound_arg, double upper_bound_arg, double fhp_arg, double shp_arg) : - shape(shape_arg), mean(mean_arg), mode(mode_arg), standard(standard_arg), lower_bound(lower_bound_arg), upper_bound(upper_bound_arg), fhp(fhp_arg), shp(shp_arg) +Prior::Prior(double mean_arg, double mode_arg, double standard_arg, double lower_bound_arg, double upper_bound_arg, double fhp_arg, double shp_arg) : + mean(mean_arg), mode(mode_arg), standard(standard_arg), lower_bound(lower_bound_arg), upper_bound(upper_bound_arg), fhp(fhp_arg), shp(shp_arg) { } diff --git a/mex/sources/estimation/Prior.hh b/mex/sources/estimation/Prior.hh index ca9ad4322..df99894b0 100644 --- a/mex/sources/estimation/Prior.hh +++ b/mex/sources/estimation/Prior.hh @@ -1,21 +1,21 @@ /* - * 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 . - */ +* 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 . +*/ /////////////////////////////////////////////////////////// // Prior.h @@ -26,39 +26,208 @@ #if !defined(Prior_8D5F562F_C831_43f3_B390_5C4EF4433756__INCLUDED_) #define Prior_8D5F562F_C831_43f3_B390_5C4EF4433756__INCLUDED_ -struct Prior -{ +#include // for beta_distribution. +#include // for gamma_distribution. +#include // for normal_distribution. +#include // for uniform_distribution. + -public: - //! probablity density functions - enum pShape +struct Prior { - Beta = 1, - Gamma = 2, - Gaussian = 3, // i.e. Normal density - Inv_gamma_1 = 4, // Inverse gamma (type 1) density - Uniform = 5, - Inv_gamma_2 = 6 //Inverse gamma (type 2) density + + public: + //! probablity density functions + enum pShape + { + Beta = 1, + Gamma = 2, + Gaussian = 3, // i.e. Normal density + Inv_gamma_1 = 4, // Inverse gamma (type 1) density + Uniform = 5, + Inv_gamma_2 = 6 //Inverse gamma (type 2) density + }; + + Prior(double mean, double mode, double standard, double lower_bound, double upper_bound, double fhp, double shp); + virtual ~Prior(); + + const double mean; + const double mode; + const double standard; + const double lower_bound; + const double upper_bound; + /** + * first shape parameter + */ + const double fhp; + /** + * second shape parameter + */ + const double shp; + virtual pShape getShape()=0; // e.g. = Beta for beta shape? + + virtual double pdf(double x) // probability density function value for x + { + std::cout << "Parent pdf undefined at parent level" << std::endl ; + return 0.0; + }; + + virtual double drand(double x) // rand for density + { + std::cout << "Parent rand undefined at parent level" << std::endl ; + return 0.0; + }; + }; - Prior(Prior::pShape shape, double mean, double mode, double standard, double lower_bound, double upper_bound, double fhp, double shp); - virtual ~Prior(); +struct BetaPrior: public Prior + { + public: + boost::math::beta_distribution distribution; - const pShape shape; - const double mean; - const double mode; - const double standard; - const double lower_bound; - const double upper_bound; - /** - * first shape parameter - */ - const double fhp; - /** - * second shape parameter - */ - const double shp; + BetaPrior( double mean, double mode, double standard, double lower_bound, double upper_bound, double fhp, double shp) + : Prior( mean, mode, standard, lower_bound, upper_bound, fhp, shp), + distribution(fhp, shp) + {}; + virtual ~BetaPrior(){}; + virtual pShape getShape(){return Prior::Beta;}; // e.g. = Beta for beta shape? -}; + virtual double pdf(double x) + { + double scalled=x; + if (lower_bound ||1.0-upper_bound) + scalled=(x- lower_bound)/(upper_bound- lower_bound); + return boost::math::pdf(distribution,scalled); + }; + + virtual double drand(double x) // rand for density + { + return 0.0; + }; + }; + +struct GammaPrior: public Prior + { + public: + boost::math::gamma_distribution distribution; + + GammaPrior( double mean, double mode, double standard, + double lower_bound, double upper_bound, double fhp, double shp) + : Prior( mean, mode, standard, lower_bound, upper_bound, fhp, shp), + distribution(fhp, shp) + {}; + virtual ~GammaPrior(){}; + virtual pShape getShape(){return Prior::Gamma;}; // e.g. = Beta for beta shape? + virtual double pdf(double x) + { + return boost::math::pdf(distribution,x- lower_bound); + }; + virtual double drand(double x) // rand for density + { + return 0.0; + }; + }; + +// X ~ IG1(s,nu) if X = sqrt(Y) where Y ~ IG2(s,nu) and Y = inv(Z) with Z ~ G(nu/2,2/s) (Gamma distribution) +// lpdfig1(x,s,n)= lpdfgam(1/(x*x),n/2,2/s)-2*log(x*x) +struct InvGamma1_Prior: public Prior + { + public: + boost::math::gamma_distribution distribution; + InvGamma1_Prior( double mean, double mode, double standard, + double lower_bound, double upper_bound, double fhp, double shp) + : Prior( mean, mode, standard, lower_bound, upper_bound, fhp, shp), + distribution(shp/2, 2/fhp) + {}; + virtual ~InvGamma1_Prior(){}; + virtual pShape getShape(){return Prior::Inv_gamma_1;}; // e.g. = Beta for beta shape? + virtual double pdf(double x) + { + double scalled=((x- lower_bound)*(x-lower_bound)); + if (x>lower_bound) + return boost::math::pdf(distribution,1/scalled) / (scalled*scalled); + else + return 0; + }; + virtual double drand(double x) // rand for density + { + return 0.0; + }; + }; + +// If x~InvGamma(a,b) , then 1/x ~Gamma(a,1/b) distribution +struct InvGamma2_Prior: public Prior + { + public: + boost::math::gamma_distribution distribution; + + InvGamma2_Prior( double mean, double mode, double standard, + double lower_bound, double upper_bound, double fhp, double shp) + : Prior( mean, mode, standard, lower_bound, upper_bound, fhp, shp), + distribution(shp/2, 2/fhp) + {}; + virtual ~InvGamma2_Prior(){}; + virtual pShape getShape(){return Prior::Inv_gamma_2;}; // e.g. = Beta for beta shape? + + virtual double pdf(double x) + { + double scalled= x - lower_bound; + if (scalled>0) + return boost::math::pdf(distribution,1/scalled) / (scalled*scalled); + else + return 0; + }; + virtual double drand(double x) // rand for density + { + return 0.0; + }; + }; + +struct GaussianPrior: public Prior + { + public: + boost::math::normal_distribution distribution; + + GaussianPrior( double mean, double mode, double standard, double lower_bound, double upper_bound, double fhp, double shp) + : Prior( mean, mode, standard, lower_bound, upper_bound, fhp, shp), + distribution(fhp, shp) //distribution(mean, standard) + {}; + virtual ~GaussianPrior(){}; + virtual pShape getShape(){return Prior::Gaussian;}; // e.g. = Beta for beta shape? + virtual double pdf(double x) + { + if (x>lower_bound && x distribution; + + UniformPrior( double mean, double mode, double standard, double lower_bound, double upper_bound, double fhp, double shp) + : Prior( mean, mode, standard, lower_bound, upper_bound, fhp, shp), + distribution(fhp, shp) //distribution(lower_bound, upper_bound) + {}; + virtual ~UniformPrior(){}; + virtual pShape getShape(){return Prior::Uniform;}; // e.g. = Beta for beta shape? + virtual double pdf(double x) + { + if (x>lower_bound && x