From b61c9e6a69c2fa82c2213b6fb1ce8053f4890132 Mon Sep 17 00:00:00 2001 From: sebastien Date: Tue, 5 Aug 2008 10:48:04 +0000 Subject: [PATCH] v4 matlab: * added preprocessor support for inverse gamma of type 2 * added support for this distribution in prior_bounds.m and rndprior.m * other cosmetic changes git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1983 ac1d8469-bf42-47a9-8791-bf33cf982152 --- matlab/draw_prior_density.m | 15 +++++++++------ matlab/prior_bounds.m | 6 ++++-- matlab/prior_draw.m | 6 ++---- matlab/rndprior.m | 17 ++++++++++++----- preprocessor/DynareBison.yy | 6 ++++-- preprocessor/DynareFlex.ll | 6 +++--- 6 files changed, 34 insertions(+), 22 deletions(-) diff --git a/matlab/draw_prior_density.m b/matlab/draw_prior_density.m index 98e5215a8..138186faf 100644 --- a/matlab/draw_prior_density.m +++ b/matlab/draw_prior_density.m @@ -44,7 +44,8 @@ p4 = bayestopt_.p4; truncprior = 10^(-3); -if pshape(indx) == 1 %/* BETA Prior */ +switch pshape(indx) + case 1 % Beta prior density = inline('((bb-x).^(b-1)).*(x-aa).^(a-1)./(beta(a,b)*(bb-aa)^(a+b-1))','x','a','b','aa','bb'); mu = (p1(indx)-p3(indx))/(p4(indx)-p3(indx)); stdd = p2(indx)/(p4(indx)-p3(indx)); @@ -57,7 +58,7 @@ if pshape(indx) == 1 %/* BETA Prior */ stepsize = (supbound-infbound)/200; abscissa = infbound:stepsize:supbound; dens = density(abscissa,a,b,aa,bb); -elseif pshape(indx) == 2 %/* GAMMA PRIOR */ + case 2 % Generalized Gamma prior mu = p1(indx)-p3(indx); b = p2(indx)^2/mu; a = mu/b; @@ -67,7 +68,7 @@ elseif pshape(indx) == 2 %/* GAMMA PRIOR */ abscissa = infbound:stepsize:supbound; dens = exp(lpdfgam(abscissa,a,b)); abscissa = abscissa + p3(indx); -elseif pshape(indx) == 3 %/* GAUSSIAN PRIOR */ + case 3 % Gaussian prior density = inline('inv(sqrt(2*pi)*b)*exp(-0.5*((x-a)/b).^2)','x','a','b'); a = p1(indx); b = p2(indx); @@ -76,7 +77,7 @@ elseif pshape(indx) == 3 %/* GAUSSIAN PRIOR */ stepsize = (supbound-infbound)/200; abscissa = infbound:stepsize:supbound; dens = density(abscissa,a,b); -elseif pshape(indx) == 4 %/* INVGAMMA PRIOR type 1 */ + case 4 % Inverse-gamma of type 1 prior density = inline('2*inv(gamma(nu/2))*(x.^(-nu-1))*((s/2)^(nu/2)).*exp(-s./(2*x.^2))','x','s','nu'); nu = p2(indx); s = p1(indx); @@ -87,7 +88,7 @@ elseif pshape(indx) == 4 %/* INVGAMMA PRIOR type 1 */ stepsize = (supbound-infbound)/200; abscissa = infbound:stepsize:supbound; dens = density(abscissa,s,nu); -elseif pshape(indx) == 5 %/* UNIFORM PRIOR */ + case 5 % Uniform prior density = inline('(x.^0)/(b-a)','x','a','b'); a = p1(indx); b = p2(indx); @@ -96,7 +97,7 @@ elseif pshape(indx) == 5 %/* UNIFORM PRIOR */ stepsize = (supbound-infbound)/200; abscissa = infbound:stepsize:supbound; dens = density(abscissa,a,b); -elseif pshape(indx) == 6 %/* INVGAMMA PRIOR type 2 */ + case 6 % Inverse-gamma of type 2 prior density = inline('inv(gamma(nu/2))*(x.^(-.5*(nu+2)))*((s/2)^(nu/2)).*exp(-s./(2*x))','x','s','nu'); nu = p2(indx); s = p1(indx); @@ -107,6 +108,8 @@ elseif pshape(indx) == 6 %/* INVGAMMA PRIOR type 2 */ stepsize = (supbound-infbound)/200; abscissa = infbound:stepsize:supbound; dens = density(abscissa,s,nu); + otherwise + error(sprintf('draw_prior_density: unknown distribution shape (index %d, type %d)', indx, pshape(indx))); end k = [1:length(dens)]; diff --git a/matlab/prior_bounds.m b/matlab/prior_bounds.m index d4ef05142..57144c523 100644 --- a/matlab/prior_bounds.m +++ b/matlab/prior_bounds.m @@ -67,8 +67,10 @@ for i=1:n case 5 bounds(i,1) = p1(i); bounds(i,2) = p2(i); + case 6 + bounds(i,1) = 1/gaminv(1-options_.prior_trunc, p2(i)/2, 2/p1(i)); + bounds(i,2) = 1/gaminv(options_.prior_trunc, p2(i)/2, 2/p1(i)); otherwise - bounds(i,1) = -Inf; - bounds(i,2) = Inf; + error(sprintf('prior_bounds: unknown distribution shape (index %d, type %d)', i, pshape(i))); end end \ No newline at end of file diff --git a/matlab/prior_draw.m b/matlab/prior_draw.m index 25b921f43..8b298bb4c 100644 --- a/matlab/prior_draw.m +++ b/matlab/prior_draw.m @@ -82,9 +82,7 @@ if init % 5: Uniform prior % p3(i) and p4(i) are used. otherwise - disp('prior_draw :: Error!') - disp('Unknown prior shape.') - return + error(sprintf('prior_draw: unknown distribution shape (index %d, type %d)', i, pshape(i))); end pdraw = zeros(npar,1); end @@ -139,6 +137,6 @@ for i = 1:npar end end otherwise - % Nothing to do here. + error(sprintf('prior_draw: unknown distribution shape (index %d, type %d)', i, pshape(i))); end end diff --git a/matlab/rndprior.m b/matlab/rndprior.m index 0ffa86b44..5b45aa8aa 100644 --- a/matlab/rndprior.m +++ b/matlab/rndprior.m @@ -39,7 +39,7 @@ for i=1:length(pmean), switch pshape(i) - case 1 %'beta' + case 1 % Beta mu = (pmean(i)-p3(i))/(p4(i)-p3(i)); stdd = p2(i)/(p4(i)-p3(i)); A = (1-mu)*mu^2/stdd^2 - mu; @@ -47,25 +47,32 @@ for i=1:length(pmean), y(1,i) = betarnd(A, B); y(1,i) = y(1,i) * (p4(i)-p3(i)) + p3(i); - case 2 %'gamma' + case 2 % Generalized gamma mu = pmean(i)-p3(i); B = p2(i)^2/mu; A = mu/B; y(1,i) = gamrnd(A, B) + p3(i); - case 3 %'normal' + case 3 % Gaussian MU = pmean(i); SIGMA = p2(i); y(1,i) = randn*SIGMA+ MU; - case 4 %'invgamma' + case 4 % Inverse-gamma type 1 nu = p2(i); s = p1(i); y(1,i) = 1/sqrt(gamrnd(nu/2, 2/s)); - case 5 %'uniform' + case 5 % Uniform y(1,i) = rand*(p2(i)-p1(i)) + p1(i); + + case 6 % Inverse-gamma type 2 + nu = p2(i); + s = p1(i); + y(1,i) = 1/gamrnd(nu/2, 2/s); + otherwise + error(sprintf('rndprior: unknown distribution shape (index %d, type %d)', i, pshape(i))); end end diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy index e2b50ad0f..85c0c0865 100644 --- a/preprocessor/DynareBison.yy +++ b/preprocessor/DynareBison.yy @@ -94,7 +94,7 @@ class ParsingDriver; %token HISTVAL HP_FILTER HP_NGRID %token INITVAL INITVAL_FILE %token INT_NUMBER -%token INV_GAMMA_PDF IRF +%token INV_GAMMA1_PDF INV_GAMMA2_PDF IRF %token KALMAN_ALGO KALMAN_TOL %token LAPLACE LCC_COMPILER LIK_ALGO LIK_INIT LINEAR LOAD_MH_FILE LOGLINEAR LU MARKOWITZ MAX %token METHOD MH_DROP MH_INIT_SCALE MH_JSCALE MH_MODE MH_NBLOCKS MH_REPLIC MH_RECOVER MIN @@ -897,10 +897,12 @@ prior : BETA_PDF { $$ = new string("2"); } | NORMAL_PDF { $$ = new string("3"); } - | INV_GAMMA_PDF + | INV_GAMMA1_PDF { $$ = new string("4"); } | UNIFORM_PDF { $$ = new string("5"); } + | INV_GAMMA2_PDF + { $$ = new string("6"); } ; value : { $$ = new string("NaN"); } diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll index f99abc416..e222fe5df 100644 --- a/preprocessor/DynareFlex.ll +++ b/preprocessor/DynareFlex.ll @@ -241,9 +241,9 @@ int sigma_e = 0; gamma_pdf {return token::GAMMA_PDF;} beta_pdf {return token::BETA_PDF;} normal_pdf {return token::NORMAL_PDF;} -inv_gamma_pdf {return token::INV_GAMMA_PDF;} -inv_gamma1_pdf {return token::INV_GAMMA_PDF;} -inv_gamma2_pdf {return token::INV_GAMMA_PDF;} +inv_gamma_pdf {return token::INV_GAMMA1_PDF;} +inv_gamma1_pdf {return token::INV_GAMMA1_PDF;} +inv_gamma2_pdf {return token::INV_GAMMA2_PDF;} uniform_pdf {return token::UNIFORM_PDF;} ; {return Dynare::parser::token_type (yytext[0]);}