From 69abb07f05cedac87a0460cb60a58fe19f0ac08f Mon Sep 17 00:00:00 2001 From: George Perendia Date: Tue, 8 Jun 2010 16:37:35 +0100 Subject: [PATCH] C++ Estimation DLL: Adding basic Uniform and Gaussian random number generators to Priors.hh and a basic test routine --- mex/sources/estimation/Prior.hh | 466 ++++++++++++++---------- mex/sources/estimation/tests/testPDF.cc | 79 ++++ 2 files changed, 346 insertions(+), 199 deletions(-) create mode 100644 mex/sources/estimation/tests/testPDF.cc diff --git a/mex/sources/estimation/Prior.hh b/mex/sources/estimation/Prior.hh index df99894b0..ba3f6a9dd 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,208 +26,276 @@ #if !defined(Prior_8D5F562F_C831_43f3_B390_5C4EF4433756__INCLUDED_) #define Prior_8D5F562F_C831_43f3_B390_5C4EF4433756__INCLUDED_ -#include // for beta_distribution. -#include // for gamma_distribution. -#include // for normal_distribution. -#include // for uniform_distribution. - +#include +#include +#include +#include + +#include // for beta_distribution. +#include // for gamma_distribution. +#include // for normal_distribution. +#include // for uniform_distribution. + +// typedef for base_uniform_generator_type +// rend48 seems better than the basic minstd_rand but +// one my later try ecuyer1988 or taus88 which are better but not yet supported in the Boost versionwe use +typedef boost::rand48 base_uniform_generator_type; struct Prior +{ + +public: + //! probablity density functions + enum pShape { - - 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; - }; - + 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 }; -struct BetaPrior: public Prior + 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 { - public: - boost::math::beta_distribution distribution; - - 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; - }; + std::cout << "Parent pdf undefined at parent level" << std::endl; + return 0.0; }; -struct GammaPrior: public Prior + virtual double + drand() // rand for density { - 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; - }; + std::cout << "Parent rand undefined at parent level" << std::endl; + 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) +}; + +struct BetaPrior : public Prior +{ +public: + boost::math::beta_distribution distribution; + + 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() // 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() // 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 +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) { - 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; - }; }; + 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() // rand for density + { + return 0.0; + }; +}; // If x~InvGamma(a,b) , then 1/x ~Gamma(a,1/b) distribution -struct InvGamma2_Prior: public Prior +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) { - public: - boost::math::gamma_distribution distribution; + }; + virtual ~InvGamma2_Prior(){}; + virtual pShape + getShape() + { + return Prior::Inv_gamma_2; + }; // e.g. = Beta for beta shape? - 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() // rand for density + { + return 0.0; + }; +}; - 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 +{ +private: + base_uniform_generator_type base_rng_type; + boost::normal_distribution rng_type; + boost::variate_generator > vrng; + +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), + rng_type(fhp, shp), // random number generator distribution type (mean, standard) + vrng(base_rng_type, rng_type), // random variate_generator + distribution(fhp, shp) //pdf 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 < upper_bound) + return boost::math::pdf(distribution, x); + else + return 0; + }; + virtual double + drand() // rand for density + { + return vrng(); + }; +}; + +struct UniformPrior : public Prior +{ +private: + base_uniform_generator_type base_rng_type; + boost::uniform_real<> rng_type; + boost::variate_generator > vrng; + +public: + boost::math::uniform_distribution 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), + rng_type(fhp, shp), // random number generator distribution type + vrng(base_rng_type, rng_type), // random variate_generator + distribution(fhp, shp) //pdf 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 < upper_bound) + return boost::math::pdf(distribution, x); + else + return 0; + }; + virtual double + drand() // rand for density + { + return vrng(); }; -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. + */ + +// test for Prior pdfs and basic random number generators + +#include "Prior.hh" + +int +main(int argc, char **argv) +{ +//BetaPrior bp( Prior::Beta, 0.1, 0.1, 0.1, 0.0, 1.0, 5.0, 1.0); + BetaPrior bp(0.1, 0.1, 0.1, 0.0, 1.0, 5.0, 1.0); + double pdf = bp.pdf(0.5); + std::cout << "beta pdf of 5,1, 0.5: " << std::setw(13) << pdf << std::endl; + BetaPrior *bpp = new BetaPrior( //Prior::Beta, + 0.1, 0.1, 0.1, 0.0, 1.0, 1.0, 5.0); + Prior *pp = bpp; + pdf = (*pp).pdf(0.1); + std::cout << "Parent (Beta) pdf of 1,5, 0.1: " << std::setw(13) << pdf << std::endl; + + GammaPrior *gpp = new GammaPrior( //Prior::Beta, + 0.1, 0.1, 0.1, 0.0, 1.0, 1.0, 5.0); + pp = gpp; + pdf = (*pp).pdf(0.1); + std::cout << "Parent (Gamma) pdf of 1,5, 0.1: " << std::setw(13) << pdf << std::endl; + + UniformPrior up(5, 5, 2, 1, 10, 20, 100); + std::cout << std::endl << "Uniform prior (5,5,2, 1,10,20,100): " << std::endl; + double ur, updf; + for (int i = 0; i < 10; i++) + { + ur = up.drand(); + updf = up.pdf(ur); + std::cout << "Uniform pdf of : " << ur << " = " << updf << std::endl; + } + + GaussianPrior gp(5, 5, 2, 1, 10, 20, 100); + std::cout << std::endl << "Gaussian prior (5,5,2, 1,10,20,100): " << std::endl; + for (int i = 0; i < 10; i++) + { + ur = gp.drand(); + updf = gp.pdf(ur); + std::cout << "Gaussian pdf of : " << ur << " = " << updf << std::endl; + } + + UniformPrior up2(5, 5, 2, 1, 100, 20, 100); + std::cout << std::endl << "Uniform prior (5,5,2, 1,100,20,100): " << std::endl; + for (int i = 0; i < 10; i++) + { + ur = up2.drand(); + updf = up2.pdf(ur); + std::cout << "Uniform pdf of : " << ur << " = " << updf << std::endl; + } + + GaussianPrior gp2(5, 5, 2, 1, 100, 5, 2); + std::cout << std::endl << "Gaussian prior (5,5,2, 1,100,5,2): " << std::endl; + for (int i = 0; i < 10; i++) + { + ur = gp2.drand(); + updf = gp2.pdf(ur); + std::cout << "Gaussian pdf of : " << ur << " = " << updf << std::endl; + } +};