From 3d9f12d769b4ad81b46b1d6d17c75ed08914db48 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Fri, 24 Jul 2015 15:08:52 +0200 Subject: [PATCH 1/5] Add check whether prior mean violates Generalized (Inverse) Gamma Performs check as for generalized beta distribution --- matlab/set_prior.m | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/matlab/set_prior.m b/matlab/set_prior.m index a72a6204d..16e27f3f3 100644 --- a/matlab/set_prior.m +++ b/matlab/set_prior.m @@ -196,6 +196,9 @@ k2 = find(isnan(bayestopt_.p4(k))); bayestopt_.p3(k(k1)) = zeros(length(k1),1); bayestopt_.p4(k(k2)) = Inf(length(k2),1); for i=1:length(k) + if (bayestopt_.p1(k(i))bayestopt_.p4(k(i))) + error(['The prior mean of ' bayestopt_.name{k(i)} ' has to be above the lower (' num2str(bayestopt_.p3(k(i))) ') bound of the Gamma prior density!']); + end if isinf(bayestopt_.p2(k(i))) error(['Infinite prior standard deviation for parameter ' bayestopt_.name{k(i)} ' is not allowed (Gamma prior)!']) end @@ -224,6 +227,9 @@ k2 = find(isnan(bayestopt_.p4(k))); bayestopt_.p3(k(k1)) = zeros(length(k1),1); bayestopt_.p4(k(k2)) = Inf(length(k2),1); for i=1:length(k) + if (bayestopt_.p1(k(i))bayestopt_.p4(k(i))) + error(['The prior mean of ' bayestopt_.name{k(i)} ' has to be above the lower (' num2str(bayestopt_.p3(k(i))) ') bound of the Inverse Gamma prior density!']); + end [bayestopt_.p6(k(i)),bayestopt_.p7(k(i))] = ... inverse_gamma_specification(bayestopt_.p1(k(i))-bayestopt_.p3(k(i)),bayestopt_.p2(k(i)),1,0) ; bayestopt_.p5(k(i)) = compute_prior_mode([ bayestopt_.p6(k(i)) , bayestopt_.p7(k(i)) , bayestopt_.p3(k(i)) ], 4) ; @@ -246,6 +252,9 @@ k2 = find(isnan(bayestopt_.p4(k))); bayestopt_.p3(k(k1)) = zeros(length(k1),1); bayestopt_.p4(k(k2)) = Inf(length(k2),1); for i=1:length(k) + if (bayestopt_.p1(k(i))bayestopt_.p4(k(i))) + error(['The prior mean of ' bayestopt_.name{k(i)} ' has to be above the lower (' num2str(bayestopt_.p3(k(i))) ') bound of the Inverse Gamma II prior density!']); + end [bayestopt_.p6(k(i)),bayestopt_.p7(k(i))] = ... inverse_gamma_specification(bayestopt_.p1(k(i))-bayestopt_.p3(k(i)),bayestopt_.p2(k(i)),2,0); bayestopt_.p5(k(i)) = compute_prior_mode([ bayestopt_.p6(k(i)) , bayestopt_.p7(k(i)) , bayestopt_.p3(k(i)) ], 6) ; From 560d9719db4114e34044a2e38a3525737faa6943 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Fri, 24 Jul 2015 15:59:40 +0200 Subject: [PATCH 2/5] Correct copy and paste mistake in prior_draw.m for inverse gamma II In case of bound violations, it was instead sampled from inverse gamma I --- matlab/prior_draw.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/prior_draw.m b/matlab/prior_draw.m index 8923be87c..9a22dc39c 100644 --- a/matlab/prior_draw.m +++ b/matlab/prior_draw.m @@ -139,7 +139,7 @@ if inverse_gamma_2_draws out_of_bound = find( (pdraw(inverse_gamma_2_index)'>ub(inverse_gamma_2_index)) | (pdraw(inverse_gamma_2_index)'ub(inverse_gamma_2_index)) | (pdraw(inverse_gamma_2_index)' Date: Fri, 24 Jul 2015 16:23:40 +0200 Subject: [PATCH 3/5] Make prior_draw.m check for uniform distribution violation of bounds Follows the logic for all other distributions --- matlab/prior_draw.m | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/matlab/prior_draw.m b/matlab/prior_draw.m index 9a22dc39c..937730748 100644 --- a/matlab/prior_draw.m +++ b/matlab/prior_draw.m @@ -1,4 +1,4 @@ -function pdraw = prior_draw(init,uniform) +function pdraw = prior_draw(init,uniform) % --*-- Unitary tests --*-- % This function generate one draw from the joint prior distribution. % % INPUTS @@ -93,6 +93,11 @@ end if uniform_draws pdraw(uniform_index) = rand(length(uniform_index),1).*(p4(uniform_index)-p3(uniform_index)) + p3(uniform_index); + out_of_bound = find( (pdraw(uniform_index)'>ub(uniform_index)) | (pdraw(uniform_index)'ub(uniform_index)) | (pdraw(uniform_index)'ub(inverse_gamma_2_index)) | (pdraw(inverse_gamma_2_index)' Date: Fri, 24 Jul 2015 16:32:05 +0200 Subject: [PATCH 4/5] Add unit tests for prior sampling --- matlab/prior_draw.m | 368 +++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 367 insertions(+), 1 deletion(-) diff --git a/matlab/prior_draw.m b/matlab/prior_draw.m index 937730748..cd0c32a50 100644 --- a/matlab/prior_draw.m +++ b/matlab/prior_draw.m @@ -15,7 +15,7 @@ function pdraw = prior_draw(init,uniform) % --*-- Unitary tests --*-- % NOTE 1. Input arguments 1 an 2 are only needed for initialization. % NOTE 2. A given draw from the joint prior distribution does not satisfy BK conditions a priori. -% Copyright (C) 2006-2010 Dynare Team +% Copyright (C) 2006-2015 Dynare Team % % This file is part of Dynare. % @@ -149,3 +149,369 @@ if inverse_gamma_2_draws end end +%@test:1 +%$ %% Initialize required structures +%$ options_.prior_trunc=0; +%$ options_.initialize_estimated_parameters_with_the_prior_mode=0; +%$ +%$ M_.dname='test'; +%$ M_.param_names = 'alp'; +%$ ndraws=100000; +%$ global estim_params_ +%$ estim_params_.var_exo = []; +%$ estim_params_.var_endo = []; +%$ estim_params_.corrx = []; +%$ estim_params_.corrn = []; +%$ estim_params_.param_vals = []; +%$ estim_params_.param_vals = [1, NaN, (-Inf), Inf, 1, 0.356, 0.02, NaN, NaN, NaN ]; +%$ +%$ %% beta +%$ estim_params_.param_vals(1,3)= -Inf;%LB +%$ estim_params_.param_vals(1,4)= +Inf;%UB +%$ estim_params_.param_vals(1,5)= 1;%Shape +%$ estim_params_.param_vals(1,6)=0.5; +%$ estim_params_.param_vals(1,7)=0.01; +%$ estim_params_.param_vals(1,8)=NaN; +%$ estim_params_.param_vals(1,9)=NaN; +%$ +%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); +%$ +%$ pdraw = prior_draw(1,0); +%$ pdraw_vec=NaN(ndraws,1); +%$ for ii=1:ndraws +%$ pdraw_vec(ii)=prior_draw(0,0); +%$ end +%$ +%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>1e-4 || abs(std(pdraw_vec)-estim_params_.param_vals(1,7))>1e-4 || any(pdraw_vec<0) || any(pdraw_vec>1) +%$ error('Beta prior wrong') +%$ end +%$ +%$ +%$ %% Gamma +%$ estim_params_.param_vals(1,3)= -Inf;%LB +%$ estim_params_.param_vals(1,4)= +Inf;%UB +%$ estim_params_.param_vals(1,5)= 2;%Shape +%$ estim_params_.param_vals(1,6)=0.5; +%$ estim_params_.param_vals(1,7)=0.01; +%$ estim_params_.param_vals(1,8)=NaN; +%$ estim_params_.param_vals(1,9)=NaN; +%$ +%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); +%$ +%$ pdraw = prior_draw(1,0); +%$ pdraw_vec=NaN(ndraws,1); +%$ for ii=1:ndraws +%$ pdraw_vec(ii)=prior_draw(0,0); +%$ end +%$ +%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>1e-4 || abs(std(pdraw_vec)-estim_params_.param_vals(1,7))>1e-4 || any(pdraw_vec<0) +%$ error('Gamma prior wrong') +%$ end +%$ +%$ %% Normal +%$ estim_params_.param_vals(1,3)= -Inf;%LB +%$ estim_params_.param_vals(1,4)= +Inf;%UB +%$ estim_params_.param_vals(1,5)= 3;%Shape +%$ estim_params_.param_vals(1,6)=0.5; +%$ estim_params_.param_vals(1,7)=0.01; +%$ estim_params_.param_vals(1,8)=NaN; +%$ estim_params_.param_vals(1,9)=NaN; +%$ +%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); +%$ +%$ pdraw = prior_draw(1,0); +%$ pdraw_vec=NaN(ndraws,1); +%$ for ii=1:ndraws +%$ pdraw_vec(ii)=prior_draw(0,0); +%$ end +%$ +%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>1e-4 || abs(std(pdraw_vec)-estim_params_.param_vals(1,7))>1e-4 +%$ error('Normal prior wrong') +%$ end +%$ +%$ %% inverse gamma distribution (type 1) +%$ estim_params_.param_vals(1,3)= -Inf;%LB +%$ estim_params_.param_vals(1,4)= +Inf;%UB +%$ estim_params_.param_vals(1,5)= 4;%Shape +%$ estim_params_.param_vals(1,6)=0.5; +%$ estim_params_.param_vals(1,7)=0.01; +%$ estim_params_.param_vals(1,8)=NaN; +%$ estim_params_.param_vals(1,9)=NaN; +%$ +%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); +%$ +%$ pdraw = prior_draw(1,0); +%$ pdraw_vec=NaN(ndraws,1); +%$ for ii=1:ndraws +%$ pdraw_vec(ii)=prior_draw(0,0); +%$ end +%$ +%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>1e-4 || abs(std(pdraw_vec)-estim_params_.param_vals(1,7))>1e-4 || any(pdraw_vec<0) +%$ error('inverse gamma distribution (type 1) prior wrong') +%$ end +%$ +%$ %% Uniform +%$ estim_params_.param_vals(1,3)= -Inf;%LB +%$ estim_params_.param_vals(1,4)= +Inf;%UB +%$ estim_params_.param_vals(1,5)= 5;%Shape +%$ estim_params_.param_vals(1,6)=0.5; +%$ estim_params_.param_vals(1,7)=sqrt(12)^(-1)*(1-0); +%$ estim_params_.param_vals(1,8)=NaN; +%$ estim_params_.param_vals(1,9)=NaN; +%$ +%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); +%$ +%$ pdraw = prior_draw(1,0); +%$ pdraw_vec=NaN(ndraws,1); +%$ for ii=1:ndraws +%$ pdraw_vec(ii)=prior_draw(0,0); +%$ end +%$ +%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>1e-2 || abs(std(pdraw_vec)-estim_params_.param_vals(1,7))>1e-3 || any(pdraw_vec<0) || any(pdraw_vec>1) +%$ error('Uniform prior wrong') +%$ end +%$ +%$ %% inverse gamma distribution (type 2) +%$ estim_params_.param_vals(1,3)= -Inf;%LB +%$ estim_params_.param_vals(1,4)= +Inf;%UB +%$ estim_params_.param_vals(1,5)= 6;%Shape +%$ estim_params_.param_vals(1,6)=0.5; +%$ estim_params_.param_vals(1,7)=0.01; +%$ estim_params_.param_vals(1,8)=NaN; +%$ estim_params_.param_vals(1,9)=NaN; +%$ +%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); +%$ +%$ pdraw = prior_draw(1,0); +%$ pdraw_vec=NaN(ndraws,1); +%$ for ii=1:ndraws +%$ pdraw_vec(ii)=prior_draw(0,0); +%$ end +%$ +%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>1e-4 || abs(std(pdraw_vec)-estim_params_.param_vals(1,7))>1e-4 || any(pdraw_vec<0) +%$ error('inverse gamma distribution (type 2) prior wrong') +%$ end +%$ +%$ +%$ %%%%%%%%%%%%%%%%%%%%%% Generalized distributions %%%%%%%%%%%%%%%%%%%%% +%$ +%$ %% beta +%$ estim_params_.param_vals(1,3)= -Inf;%LB +%$ estim_params_.param_vals(1,4)= +Inf;%UB +%$ estim_params_.param_vals(1,5)= 1;%Shape +%$ estim_params_.param_vals(1,6)=1.5; +%$ estim_params_.param_vals(1,7)=0.01; +%$ estim_params_.param_vals(1,8)=1; +%$ estim_params_.param_vals(1,9)=3; +%$ +%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); +%$ +%$ pdraw = prior_draw(1,0); +%$ pdraw_vec=NaN(ndraws,1); +%$ for ii=1:ndraws +%$ pdraw_vec(ii)=prior_draw(0,0); +%$ end +%$ +%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>1e-4 || abs(std(pdraw_vec)-estim_params_.param_vals(1,7))>1e-4 || any(pdraw_vecestim_params_.param_vals(1,4)) +%$ error('Beta prior wrong') +%$ end +%$ +%$ %% Gamma +%$ estim_params_.param_vals(1,3)= -Inf;%LB +%$ estim_params_.param_vals(1,4)= +Inf;%UB +%$ estim_params_.param_vals(1,5)= 2;%Shape +%$ estim_params_.param_vals(1,6)=1.5; +%$ estim_params_.param_vals(1,7)=0.01; +%$ estim_params_.param_vals(1,8)=1; +%$ estim_params_.param_vals(1,9)=NaN; +%$ +%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); +%$ +%$ pdraw = prior_draw(1,0); +%$ pdraw_vec=NaN(ndraws,1); +%$ for ii=1:ndraws +%$ pdraw_vec(ii)=prior_draw(0,0); +%$ end +%$ +%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>1e-4 || abs(std(pdraw_vec)-estim_params_.param_vals(1,7))>1e-4 || any(pdraw_vec1e-4 || abs(std(pdraw_vec)-estim_params_.param_vals(1,7))>1e-4 || any(pdraw_vec1e-4 || abs(std(pdraw_vec)-estim_params_.param_vals(1,7))>1e-4 || any(pdraw_vecestim_params_.param_vals(1,4)) +%$ error('Uniform prior wrong') +%$ end +%$ +%$ %% inverse gamma distribution (type 2) +%$ estim_params_.param_vals(1,3)= -Inf;%LB +%$ estim_params_.param_vals(1,4)= +Inf;%UB +%$ estim_params_.param_vals(1,5)= 6;%Shape +%$ estim_params_.param_vals(1,6)=1.5; +%$ estim_params_.param_vals(1,7)=0.01; +%$ estim_params_.param_vals(1,8)=1; +%$ estim_params_.param_vals(1,9)=NaN; +%$ +%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); +%$ +%$ pdraw = prior_draw(1,0); +%$ pdraw_vec=NaN(ndraws,1); +%$ for ii=1:ndraws +%$ pdraw_vec(ii)=prior_draw(0,0); +%$ end +%$ +%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>1e-4 || abs(std(pdraw_vec)-estim_params_.param_vals(1,7))>1e-4 || any(pdraw_vec5e-3 || any(pdraw_vecbounds.ub) +%$ error('Beta prior wrong') +%$ end +%$ +%$ %% Gamma +%$ estim_params_.param_vals(1,3)= -Inf;%LB +%$ estim_params_.param_vals(1,4)= +Inf;%UB +%$ estim_params_.param_vals(1,5)= 2;%Shape +%$ estim_params_.param_vals(1,6)=1.5; +%$ estim_params_.param_vals(1,7)=0.01; +%$ estim_params_.param_vals(1,8)=1; +%$ estim_params_.param_vals(1,9)=NaN; +%$ +%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); +%$ bounds = prior_bounds(bayestopt_,options_)'; +%$ +%$ pdraw = prior_draw(1,0); +%$ pdraw_vec=NaN(ndraws,1); +%$ for ii=1:ndraws +%$ pdraw_vec(ii)=prior_draw(0,0); +%$ end +%$ +%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>5e-3 || any(pdraw_vecbounds.ub) +%$ error('Gamma prior wrong') +%$ end +%$ +%$ %% inverse gamma distribution (type 1) +%$ estim_params_.param_vals(1,3)= -Inf;%LB +%$ estim_params_.param_vals(1,4)= +Inf;%UB +%$ estim_params_.param_vals(1,5)= 4;%Shape +%$ estim_params_.param_vals(1,6)=1.5; +%$ estim_params_.param_vals(1,7)=0.01; +%$ estim_params_.param_vals(1,8)=1; +%$ estim_params_.param_vals(1,9)=NaN; +%$ +%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); +%$ bounds = prior_bounds(bayestopt_,options_)'; +%$ +%$ pdraw = prior_draw(1,0); +%$ pdraw_vec=NaN(ndraws,1); +%$ for ii=1:ndraws +%$ pdraw_vec(ii)=prior_draw(0,0); +%$ end +%$ +%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>5e-3 || any(pdraw_vecbounds.ub) +%$ error('inverse gamma distribution (type 1) prior wrong') +%$ end +%$ +%$ %% Uniform +%$ estim_params_.param_vals(1,3)= -Inf;%LB +%$ estim_params_.param_vals(1,4)= +Inf;%UB +%$ estim_params_.param_vals(1,5)= 5;%Shape +%$ estim_params_.param_vals(1,6)=1.5; +%$ estim_params_.param_vals(1,7)=0.01; +%$ estim_params_.param_vals(1,8)=NaN; +%$ estim_params_.param_vals(1,9)=NaN; +%$ +%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); +%$ bounds = prior_bounds(bayestopt_,options_)'; +%$ +%$ pdraw = prior_draw(1,0); +%$ pdraw_vec=NaN(ndraws,1); +%$ for ii=1:ndraws +%$ pdraw_vec(ii)=prior_draw(0,0); +%$ end +%$ +%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>5e-3 || any(pdraw_vecbounds.ub) +%$ error('Uniform prior wrong') +%$ end +%$ +%$ +%$ %% inverse gamma distribution (type 2) +%$ estim_params_.param_vals(1,3)= -Inf;%LB +%$ estim_params_.param_vals(1,4)= +Inf;%UB +%$ estim_params_.param_vals(1,5)= 6;%Shape +%$ estim_params_.param_vals(1,6)=1.5; +%$ estim_params_.param_vals(1,7)=0.01; +%$ estim_params_.param_vals(1,8)=1; +%$ estim_params_.param_vals(1,9)=NaN; +%$ +%$ [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_); +%$ bounds = prior_bounds(bayestopt_,options_)'; +%$ +%$ pdraw = prior_draw(1,0); +%$ pdraw_vec=NaN(ndraws,1); +%$ for ii=1:ndraws +%$ pdraw_vec(ii)=prior_draw(0,0); +%$ end +%$ +%$ if abs(mean(pdraw_vec)-estim_params_.param_vals(1,6))>5e-3 || any(pdraw_vecbounds.ub) +%$ error('inverse gamma distribution (type 2) prior wrong') +%$ end +%$ +%@eof:1 From 19fdb309ac18befd9932599047f823a126349fe3 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Mon, 27 Jul 2015 10:19:02 +0200 Subject: [PATCH 5/5] Fix computation of prior mode for uniform distribution While this code is never actually used, we return the prior mean, which is directly given by the first hyperparameter. --- matlab/distributions/compute_prior_mode.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/matlab/distributions/compute_prior_mode.m b/matlab/distributions/compute_prior_mode.m index 50df6a168..e6d35920c 100644 --- a/matlab/distributions/compute_prior_mode.m +++ b/matlab/distributions/compute_prior_mode.m @@ -22,7 +22,7 @@ function m = compute_prior_mode(hyperparameters,shape) % [3] The uniform distribution has an infinity of modes. In this case the function returns the prior mean. % [4] For the beta distribution we can have 1, 2 or an infinity of modes. -% Copyright (C) 2009 Dynare Team +% Copyright (C) 2009-2015 Dynare Team % % This file is part of Dynare. % @@ -74,7 +74,7 @@ switch shape m = m + hyperparameters(3); end case 5 - m = .5*(hyperparameters(2)-hyperparameters(1)) ; + m = hyperparameters(1); case 6 % s = hyperparameters(1) % nu = hyperparameters(2)