From 015513380f7735ed665f3acd56712d4b86b80006 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Ry=C3=BBk=29?= Date: Wed, 26 Apr 2023 10:34:25 +0200 Subject: [PATCH] Add new class for priors. This commit only changes the routine used to draw random deviates from the prior distribution, without relying on persistent variables (which allows parallelisation). --- matlab/@dprior/dprior.m | 376 ++++++++++++++++++++++ matlab/GetOneDraw.m | 21 +- matlab/PosteriorIRF_core1.m | 13 +- matlab/dynare_identification.m | 21 +- matlab/execute_prior_posterior_function.m | 56 ++-- matlab/isprior.m | 22 ++ matlab/particles | 2 +- matlab/posterior_sampler_core.m | 4 +- matlab/posterior_sampler_initialization.m | 13 +- matlab/prior_posterior_statistics_core.m | 17 +- matlab/prior_sampler.m | 14 +- matlab/slice_sampler.m | 8 +- 12 files changed, 485 insertions(+), 82 deletions(-) create mode 100644 matlab/@dprior/dprior.m create mode 100644 matlab/isprior.m diff --git a/matlab/@dprior/dprior.m b/matlab/@dprior/dprior.m new file mode 100644 index 000000000..94c38e9e8 --- /dev/null +++ b/matlab/@dprior/dprior.m @@ -0,0 +1,376 @@ +classdef dprior + + properties + p6 = []; % Prior first hyperparameter. + p7 = []; % Prior second hyperparameter. + p3 = []; % Lower bound of the prior support. + p4 = []; % Upper bound of the prior support. + lb = []; % Truncated prior lower bound. + ub = []; % Truncated prior upper bound. + uniform_index = []; % Index for the uniform priors. + gaussian_index = []; % Index for the gaussian priors. + gamma_index = []; % Index for the gamma priors. + beta_index = []; % Index for the beta priors. + inverse_gamma_1_index = []; % Index for the inverse gamma type 1 priors. + inverse_gamma_2_index = []; % Index for the inverse gamma type 2 priors. + weibull_index = []; % Index for the weibull priors. + uniform_draws = false; + gaussian_draws = false; + gamma_draws = false; + beta_draws = false; + inverse_gamma_1_draws = false; + inverse_gamma_2_draws = false; + weibull_draws = false; + end + + methods + + function o = dprior(BayesInfo, PriorTrunc, Uniform) + % Class constructor. + % + % INPUTS + % - BayesInfo [struct] Informations about the prior distribution, aka bayestopt_. + % - PriorTrunc [double] scalar, probability mass to be excluded, aka options_.prior_trunc + % - Uniform [logical] scalar, produce uniform random deviates on the prior support. + % + % OUTPUTS + % - o [dprior] scalar, prior object. + % + % REQUIREMENTS + % None. + o.p6 = BayesInfo.p6; + o.p7 = BayesInfo.p7; + o.p3 = BayesInfo.p3; + o.p4 = BayesInfo.p4; + bounds = prior_bounds(BayesInfo, PriorTrunc); + o.lb = bounds.lb; + o.ub = bounds.ub; + if nargin>2 && Uniform + prior_shape = repmat(5, length(o.p6), 1); + else + prior_shape = BayesInfo.pshape; + end + o.beta_index = find(prior_shape==1); + if ~isempty(o.beta_index) + o.beta_draws = true; + end + o.gamma_index = find(prior_shape==2); + if ~isempty(o.gamma_index) + o.gamma_draws = true; + end + o.gaussian_index = find(prior_shape==3); + if ~isempty(o.gaussian_index) + o.gaussian_draws = true; + end + o.inverse_gamma_1_index = find(prior_shape==4); + if ~isempty(o.inverse_gamma_1_index) + o.inverse_gamma_1_draws = true; + end + o.uniform_index = find(prior_shape==5); + if ~isempty(o.uniform_index) + o.uniform_draws = true; + end + o.inverse_gamma_2_index = find(prior_shape==6); + if ~isempty(o.inverse_gamma_2_index) + o.inverse_gamma_2_draws = true; + end + o.weibull_index = find(prior_shape==8); + if ~isempty(o.weibull_index) + o.weibull_draws = true; + end + end + + function p = draw(o) + % Return a random draw from the prior distribution. + % + % INPUTS + % - o [dprior] + % + % OUTPUTS + % - p [double] m×1 vector, random draw from the prior distribution (m is the number of estimated parameters). + % + % REMARKS + % None. + % + % EXAMPLE + % + % >> Prior = dprior(bayestopt_, options_.prior_trunc); + % >> d = Prior.draw() + p = NaN(rows(o.lb), 1); + if o.uniform_draws + p(o.uniform_index) = rand(length(o.uniform_index),1).*(o.p4(o.uniform_index)-o.p3(o.uniform_index)) + o.p3(o.uniform_index); + out_of_bound = find( (p(o.uniform_index)>o.ub(o.uniform_index)) | (p(o.uniform_index)o.ub(o.uniform_index)) | (p(o.uniform_index)o.ub(o.gaussian_index)) | (p(o.gaussian_index)o.ub(o.gaussian_index)) | (p(o.gaussian_index)o.ub(o.gamma_index)) | (p(o.gamma_index)o.ub(o.gamma_index)) | (p(o.gamma_index)o.ub(o.beta_index)) | (p(o.beta_index)o.ub(o.beta_index)) | (p(o.beta_index)o.ub(o.inverse_gamma_1_index)) | (p(o.inverse_gamma_1_index)o.ub(inverse_gamma_1_index)) | (p(o.inverse_gamma_1_index)o.ub(o.inverse_gamma_2_index)) | (p(o.inverse_gamma_2_index)o.ub(o.inverse_gamma_2_index)) | (p(o.inverse_gamma_2_index)o.ub(o.weibull_index)) | (p(o.weibull_index)o.ub(o.weibull_index)) | (p(o.weibull_index)> Prior = dprior(bayestopt_, options_.prior_trunc); + % >> Prior.draws(1e6) + P = NaN(rows(o.lb), 1); + parfor i=1:n + P(:,i) = draw(o); + end + end + + end % methods +end % classdef --*-- Unit tests --*-- + +%@test:1 +%$ % Fill global structures with required fields... +%$ prior_trunc = 1e-10; +%$ p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1); % Prior shape +%$ p1 = .4*ones(14,1); % Prior mean +%$ p2 = .2*ones(14,1); % Prior std. +%$ p3 = NaN(14,1); +%$ p4 = NaN(14,1); +%$ p5 = NaN(14,1); +%$ p6 = NaN(14,1); +%$ p7 = NaN(14,1); +%$ +%$ for i=1:14 +%$ switch p0(i) +%$ case 1 +%$ % Beta distribution +%$ p3(i) = 0; +%$ p4(i) = 1; +%$ [p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i)); +%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 1); +%$ case 2 +%$ % Gamma distribution +%$ p3(i) = 0; +%$ p4(i) = Inf; +%$ [p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i)); +%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 2); +%$ case 3 +%$ % Normal distribution +%$ p3(i) = -Inf; +%$ p4(i) = Inf; +%$ p6(i) = p1(i); +%$ p7(i) = p2(i); +%$ p5(i) = p1(i); +%$ case 4 +%$ % Inverse Gamma (type I) distribution +%$ p3(i) = 0; +%$ p4(i) = Inf; +%$ [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false); +%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 4); +%$ case 5 +%$ % Uniform distribution +%$ [p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i)); +%$ p3(i) = p6(i); +%$ p4(i) = p7(i); +%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 5); +%$ case 6 +%$ % Inverse Gamma (type II) distribution +%$ p3(i) = 0; +%$ p4(i) = Inf; +%$ [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false); +%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 6); +%$ case 8 +%$ % Weibull distribution +%$ p3(i) = 0; +%$ p4(i) = Inf; +%$ [p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i)); +%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 8); +%$ otherwise +%$ error('This density is not implemented!') +%$ end +%$ end +%$ +%$ BayesInfo.pshape = p0; +%$ BayesInfo.p1 = p1; +%$ BayesInfo.p2 = p2; +%$ BayesInfo.p3 = p3; +%$ BayesInfo.p4 = p4; +%$ BayesInfo.p5 = p5; +%$ BayesInfo.p6 = p6; +%$ BayesInfo.p7 = p7; +%$ +%$ ndraws = 1e5; +%$ m0 = BayesInfo.p1; %zeros(14,1); +%$ v0 = diag(BayesInfo.p2.^2); %zeros(14); +%$ +%$ % Call the tested routine +%$ try +%$ % Instantiate dprior object +%$ o = dprior(BayesInfo, prior_trunc, false); +%$ % Do simulations in a loop and estimate recursively the mean and the variance. +%$ for i = 1:ndraws +%$ draw = o.draw(); +%$ m1 = m0 + (draw-m0)/i; +%$ m2 = m1*m1'; +%$ v0 = v0 + ((draw*draw'-m2-v0) + (i-1)*(m0*m0'-m2'))/i; +%$ m0 = m1; +%$ end +%$ t(1) = true; +%$ catch +%$ t(1) = false; +%$ end +%$ +%$ if t(1) +%$ t(2) = all(abs(m0-BayesInfo.p1)<3e-3); +%$ t(3) = all(all(abs(v0-diag(BayesInfo.p2.^2))<5e-3)); +%$ end +%$ T = all(t); +%@eof:1 + +%@test:2 +%$ % Fill global structures with required fields... +%$ prior_trunc = 1e-10; +%$ p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1); % Prior shape +%$ p1 = .4*ones(14,1); % Prior mean +%$ p2 = .2*ones(14,1); % Prior std. +%$ p3 = NaN(14,1); +%$ p4 = NaN(14,1); +%$ p5 = NaN(14,1); +%$ p6 = NaN(14,1); +%$ p7 = NaN(14,1); +%$ +%$ for i=1:14 +%$ switch p0(i) +%$ case 1 +%$ % Beta distribution +%$ p3(i) = 0; +%$ p4(i) = 1; +%$ [p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i)); +%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 1); +%$ case 2 +%$ % Gamma distribution +%$ p3(i) = 0; +%$ p4(i) = Inf; +%$ [p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i)); +%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 2); +%$ case 3 +%$ % Normal distribution +%$ p3(i) = -Inf; +%$ p4(i) = Inf; +%$ p6(i) = p1(i); +%$ p7(i) = p2(i); +%$ p5(i) = p1(i); +%$ case 4 +%$ % Inverse Gamma (type I) distribution +%$ p3(i) = 0; +%$ p4(i) = Inf; +%$ [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false); +%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 4); +%$ case 5 +%$ % Uniform distribution +%$ [p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i)); +%$ p3(i) = p6(i); +%$ p4(i) = p7(i); +%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 5); +%$ case 6 +%$ % Inverse Gamma (type II) distribution +%$ p3(i) = 0; +%$ p4(i) = Inf; +%$ [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false); +%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 6); +%$ case 8 +%$ % Weibull distribution +%$ p3(i) = 0; +%$ p4(i) = Inf; +%$ [p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i)); +%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 8); +%$ otherwise +%$ error('This density is not implemented!') +%$ end +%$ end +%$ +%$ BayesInfo.pshape = p0; +%$ BayesInfo.p1 = p1; +%$ BayesInfo.p2 = p2; +%$ BayesInfo.p3 = p3; +%$ BayesInfo.p4 = p4; +%$ BayesInfo.p5 = p5; +%$ BayesInfo.p6 = p6; +%$ BayesInfo.p7 = p7; +%$ +%$ ndraws = 1e5; +%$ +%$ % Call the tested routine +%$ try +%$ % Instantiate dprior object. +%$ o = dprior(BayesInfo, prior_trunc, false); +%$ X = o.draws(ndraws); +%$ m = mean(X, 2); +%$ v = var(X, 0, 2); +%$ t(1) = true; +%$ catch +%$ t(1) = false; +%$ end +%$ +%$ if t(1) +%$ t(2) = all(abs(m-BayesInfo.p1)<3e-3); +%$ t(3) = all(all(abs(v-BayesInfo.p2.^2)<5e-3)); +%$ end +%$ T = all(t); +%@eof:2 diff --git a/matlab/GetOneDraw.m b/matlab/GetOneDraw.m index ffd5665c9..8c92a84ba 100644 --- a/matlab/GetOneDraw.m +++ b/matlab/GetOneDraw.m @@ -1,5 +1,5 @@ -function [xparams, logpost] = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_) -% function [xparams, logpost] = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_) +function [xparams, logpost] = GetOneDraw(distribution, M_, estim_params_, oo_, options_, bayestopt_) + % draws one parameter vector and its posterior from MCMC or the prior % % INPUTS @@ -18,7 +18,7 @@ function [xparams, logpost] = GetOneDraw(type,M_,estim_params_,oo_,options_,baye % SPECIAL REQUIREMENTS % none -% Copyright © 2005-2017 Dynare Team +% Copyright © 2005-2023 Dynare Team % % This file is part of Dynare. % @@ -35,12 +35,13 @@ function [xparams, logpost] = GetOneDraw(type,M_,estim_params_,oo_,options_,baye % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -switch type - case 'posterior' - [xparams, logpost] = metropolis_draw(0); - case 'prior' - xparams = prior_draw(); +if isprior(distribution) + xparams = distribution.draw(); if nargout>1 - logpost = evaluate_posterior_kernel(xparams',M_,estim_params_,oo_,options_,bayestopt_); + logpost = evaluate_posterior_kernel(xparams, M_, estim_params_, oo_, options_, bayestopt_); end -end \ No newline at end of file +elseif ischar(distribution) && strcmpi(distribution, 'posterior') + [xparams, logpost] = metropolis_draw(0); +else + error('GetOneDraw:: Wrong inputs.') +end diff --git a/matlab/PosteriorIRF_core1.m b/matlab/PosteriorIRF_core1.m index 211870470..9b5122d14 100644 --- a/matlab/PosteriorIRF_core1.m +++ b/matlab/PosteriorIRF_core1.m @@ -23,7 +23,7 @@ function myoutput=PosteriorIRF_core1(myinputs,fpar,B,whoiam, ThisMatlab) % SPECIAL REQUIREMENTS. % None. % -% Copyright © 2006-2019 Dynare Team +% Copyright © 2006-2023 Dynare Team % % This file is part of Dynare. % @@ -134,12 +134,17 @@ end % Parallel 'while' very good!!! stock_param=zeros(MAX_nruns,npar); stock_irf_dsge=zeros(options_.irf,nvar,M_.exo_nbr,MAX_nirfs_dsge); + +if strcmp(type, 'prior') + Prior = dprior(bayestopt_, options_.prior_trunc); +end + while fpar 0) if options_ident.prior_range %sample uniformly from prior ranges (overwrite prior specification) - prior_draw(bayestopt_, options_.prior_trunc, true); + Prior = dprior(bayestopt_, options_.prior_trunc, true); else %sample from prior distributions - prior_draw(bayestopt_, options_.prior_trunc, false); + Prior = dprior(bayestopt_, options_.prior_trunc, false); end else options_ident.prior_mc = 1; %only one single point @@ -495,7 +496,7 @@ if iload <=0 kk=0; while kk<50 && info(1) kk=kk+1; - params = prior_draw(); + params = Prior.draw(); options_ident.tittxt = 'Random_prior_params'; %title text for graphs and figures % perform identification analysis [ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, derivatives_info_point, info, error_indicator_point] = ... @@ -548,7 +549,7 @@ if iload <=0 if external_sample params = pdraws0(iteration+1,:); % loaded draws else - params = prior_draw(); % new random draw from prior + params = Prior.draw(); % new random draw from prior end options_ident.tittxt = []; % clear title text for graphs and figures % run identification analysis @@ -692,7 +693,7 @@ if iload <=0 end % store results for minimal system - if ~options_MC.no_identification_minimal + if ~options_MC.no_identification_minimal if ~error_indicator.identification_minimal STO_dMINIMAL(:,:,run_index) = ide_minimal.dMINIMAL; IDE_MINIMAL.cond(iteration,1) = ide_minimal.cond; @@ -901,7 +902,7 @@ if SampleSize > 1 store_nodisplay = options_.nodisplay; options_.nodisplay = 1; % HIGHEST CONDITION NUMBER - [~, jmax] = max(IDE_MOMENTS.cond); + [~, jmax] = max(IDE_MOMENTS.cond); tittxt = 'Draw with HIGHEST condition number'; fprintf('\nTesting %s.\n',tittxt); if ~iload diff --git a/matlab/execute_prior_posterior_function.m b/matlab/execute_prior_posterior_function.m index 00cbe64e7..3ac13f177 100644 --- a/matlab/execute_prior_posterior_function.m +++ b/matlab/execute_prior_posterior_function.m @@ -17,7 +17,7 @@ function oo_=execute_prior_posterior_function(posterior_function_name,M_,options % OUTPUTS % oo_ [structure] Matlab/Octave structure gathering the results (initialized by dynare, see @ref{oo_}). -% Copyright © 2013-2015 Dynare Team +% Copyright © 2013-2023 Dynare Team % % This file is part of Dynare. % @@ -47,21 +47,21 @@ end %Create function handle functionhandle=str2func(posterior_function_name); -prior = true; n_draws=options_.sampling_draws; -% Get informations about the _posterior_draws files. + if strcmpi(type,'posterior') - %% discard first mh_drop percent of the draws: + % Get informations about the _posterior_draws files. + % discard first mh_drop percent of the draws: CutSample(M_, options_, estim_params_); - %% initialize metropolis draws - options_.sub_draws=n_draws; %set draws for sampling; changed value is not returned to base workspace - [error_flag,~,options_]= metropolis_draw(1,options_,estim_params_,M_); + % initialize metropolis draws + options_.sub_draws = n_draws; % set draws for sampling; changed value is not returned to base workspace + [error_flag, ~, options_] = metropolis_draw(1, options_, estim_params_, M_); if error_flag error('EXECUTE_POSTERIOR_FUNCTION: The draws could not be initialized') end - n_draws=options_.sub_draws; - prior = false; + n_draws = options_.sub_draws; elseif strcmpi(type,'prior') + % Get informations about the prior distribution. if isempty(bayestopt_) if ~isempty(estim_params_) && ~(isfield(estim_params_,'nvx') && (size(estim_params_.var_exo,1)+size(estim_params_.var_endo,1)+size(estim_params_.corrx,1)+size(estim_params_.corrn,1)+size(estim_params_.param_vals,1))==0) [xparam1,estim_params_,bayestopt_,lb,ub,M_] = set_prior(estim_params_,M_,options_); @@ -72,38 +72,36 @@ elseif strcmpi(type,'prior') if exist([M_.fname '_prior_restrictions.m']) warning('prior_function currently does not support endogenous prior restrictions. They will be ignored. Consider using a prior_function with nobs=1.') end - prior_draw(bayestopt_, options_.prior_trunc); + Prior = dprior(bayestopt_, options_.prior_trunc); else error('EXECUTE_POSTERIOR_FUNCTION: Unknown type!') end -%get draws for later use -first_draw=GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_); -parameter_mat=NaN(n_draws,length(first_draw)); -parameter_mat(1,:)=first_draw; -for draw_iter=2:n_draws - parameter_mat(draw_iter,:) = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_); +if strcmpi(type, 'prior') + parameter_mat = Prior.draws(n_draws); +else + parameter_mat = NaN(length(bayestopt_.p6), n_draws); + for i = 1:n_draws + parameter_mat(:,i) = GetOneDraw(type, M_, estim_params_, oo_, options_, bayestopt_); + end end -% get output size +% Get output size try - junk=functionhandle(parameter_mat(1,:),M_,options_,oo_,estim_params_,bayestopt_,dataset_,dataset_info); + junk = functionhandle(parameter_mat(:,1), M_, options_, oo_, estim_params_, bayestopt_, dataset_, dataset_info); catch err fprintf('\nEXECUTE_POSTERIOR_FUNCTION: Execution of prior/posterior function led to an error. Execution cancelled.\n') rethrow(err) end -%initialize cell with number of columns -results_cell=cell(n_draws,size(junk,2)); +% Initialize cell with number of columns +results_cell = cell(n_draws, columns(junk)); -%% compute function on draws -for draw_iter = 1:n_draws - M_ = set_all_parameters(parameter_mat(draw_iter,:),estim_params_,M_); - [results_cell(draw_iter,:)]=functionhandle(parameter_mat(draw_iter,:),M_,options_,oo_,estim_params_,bayestopt_,dataset_,dataset_info); +% Evaluate function on each draw +for i = 1:n_draws + M_ = set_all_parameters(parameter_mat(:,i), estim_params_, M_); + [results_cell(i,:)] = functionhandle(parameter_mat(:,i), M_, options_, oo_, estim_params_, bayestopt_, dataset_, dataset_info); end -if prior - oo_.prior_function_results = results_cell; -else - oo_.posterior_function_results = results_cell; -end +% Save results under oo_ +oo_.(sprintf('%s_function_results', type)) = results_cell; diff --git a/matlab/isprior.m b/matlab/isprior.m new file mode 100644 index 000000000..a7a014b14 --- /dev/null +++ b/matlab/isprior.m @@ -0,0 +1,22 @@ +function b = isprior(o) + +% Return True iff o is a dprior object. + +% Copyright © 2023 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 . + +b = isa(o, 'dprior'); diff --git a/matlab/particles b/matlab/particles index 7926e75aa..bb16df7f2 160000 --- a/matlab/particles +++ b/matlab/particles @@ -1 +1 @@ -Subproject commit 7926e75aa525b73c91bbb944d19556c43a8fbbe6 +Subproject commit bb16df7f2e6ddbc81e23b9661ee100f5e8dc675a diff --git a/matlab/posterior_sampler_core.m b/matlab/posterior_sampler_core.m index 7925f482b..0722858c4 100644 --- a/matlab/posterior_sampler_core.m +++ b/matlab/posterior_sampler_core.m @@ -36,7 +36,7 @@ function myoutput = posterior_sampler_core(myinputs,fblck,nblck,whoiam, ThisMatl % See the comments in the posterior_sampler.m funtion. -% Copyright © 2006-2017 Dynare Team +% Copyright © 2006-2023 Dynare Team % % This file is part of Dynare. % @@ -86,8 +86,6 @@ oo_ = myinputs.oo_; if whoiam % initialize persistent variables in priordens() priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, bayestopt_.p3,bayestopt_.p4,1); - % initialize persistent variables in prior_draw() - prior_draw(bayestopt_,options_.prior_trunc); end MetropolisFolder = CheckPath('metropolis',M_.dname); diff --git a/matlab/posterior_sampler_initialization.m b/matlab/posterior_sampler_initialization.m index 52ba42060..60a909d2c 100644 --- a/matlab/posterior_sampler_initialization.m +++ b/matlab/posterior_sampler_initialization.m @@ -38,7 +38,7 @@ function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npa % SPECIAL REQUIREMENTS % None. -% Copyright © 2006-2017 Dynare Team +% Copyright © 2006-2023 Dynare Team % % This file is part of Dynare. % @@ -113,7 +113,7 @@ if ~options_.load_mh_file && ~options_.mh_recover fprintf(fidlog,[' Number of simulations per block: ' int2str(nruns(1)) '\n']); fprintf(fidlog,' \n'); if isempty(d) - prior_draw(bayestopt_,options_.prior_trunc); + Prior = dprior(bayestopt_, options_.prior_trunc); end if options_.mh_initialize_from_previous_mcmc.status PreviousFolder0 = options_.mh_initialize_from_previous_mcmc.directory; @@ -132,7 +132,7 @@ if ~options_.load_mh_file && ~options_.mh_recover MetropolisFolder0 = fileparts(RecordFile0); PreviousFolder0=fileparts(MetropolisFolder0); [~, ModelName0]=fileparts(PreviousFolder0); - else + else %% check for proper filesep char in user defined paths PreviousFolder0=strrep(PreviousFolder0,'\',filesep); MetropolisFolder0 = [PreviousFolder0 filesep 'metropolis']; @@ -195,7 +195,7 @@ if ~options_.load_mh_file && ~options_.mh_recover trial = 1; while validate == 0 && trial <= 10 if isempty(d) - candidate = prior_draw(); + candidate = Prior.draw(); else if isfield(options_,'mh_init_scale') if trial==1 @@ -228,7 +228,7 @@ if ~options_.load_mh_file && ~options_.mh_recover if options_.nointeractive disp(['Estimation::mcmc: I reduce mh_init_scale_factor by 10 percent:']) if isfield(options_,'mh_init_scale') - options_.mh_init_scale = .9*options_.mh_init_scale; + options_.mh_init_scale = .9*options_.mh_init_scale; fprintf('Estimation::mcmc: Parameter mh_init_scale is now equal to %f.\n',options_.mh_init_scale) else options_.mh_init_scale_factor = .9*options_.mh_init_scale_factor; @@ -301,7 +301,6 @@ if ~options_.load_mh_file && ~options_.mh_recover if options_.mh_initialize_from_previous_mcmc.status record.InitialSeeds = record0.LastSeeds; else - for j=1:NumberOfBlocks % we set a different seed for the random generator for each block then we record the corresponding random generator state (vector) set_dynare_seed(options_.DynareRandomStreams.seed+j); @@ -502,7 +501,7 @@ elseif options_.mh_recover % How many mh-files are saved in this block? ExistingDrawsInLastMCFile=zeros(NumberOfBlocks,1); %initialize: no MCMC draws of current MCMC are in file from last run % Check whether last present file is a file included in the last MCMC run - + update_record=0; for k=1:NumberOfBlocks FirstBlock = k; diff --git a/matlab/prior_posterior_statistics_core.m b/matlab/prior_posterior_statistics_core.m index 7d04d2845..64cf7685a 100644 --- a/matlab/prior_posterior_statistics_core.m +++ b/matlab/prior_posterior_statistics_core.m @@ -30,7 +30,7 @@ function myoutput=prior_posterior_statistics_core(myinputs,fpar,B,whoiam, ThisMa % SPECIAL REQUIREMENTS. % None. -% Copyright © 2005-2020 Dynare Team +% Copyright © 2005-2023 Dynare Team % % This file is part of Dynare. % @@ -109,8 +109,8 @@ if ~strcmpi(type,'prior') logpost=myinputs.logpost; end else - x=(prior_draw(bayestopt_,options_.prior_trunc))'; - options_=select_qz_criterium_value(options_); + Prior = dprior(bayestopt_,options_.prior_trunc); + options_ = select_qz_criterium_value(options_); end if whoiam Parallel=myinputs.Parallel; @@ -180,9 +180,9 @@ end if naK stock_filter_step_ahead =NaN(length(options_.filter_step_ahead),endo_nbr,gend+max(options_.filter_step_ahead),MAX_naK); end -stock_param = NaN(MAX_nruns,size(x,2)); -stock_logpo = NaN(MAX_nruns,1); -stock_ys = NaN(MAX_nruns,endo_nbr); +stock_param = NaN(MAX_nruns, length(bayestopt_.p6)); +stock_logpo = NaN(MAX_nruns, 1); +stock_ys = NaN(MAX_nruns, endo_nbr); if filter_covariance stock_filter_covariance = zeros(endo_nbr,endo_nbr,gend+1,MAX_filter_covariance); end @@ -196,7 +196,8 @@ for b=fpar:B iter=1; logpo=[]; while (isempty(logpo) || isinf(logpo)) && iter<1000 - [deep, logpo] = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_); + deep = Prior.draw(); + logpo = evaluate_posterior_kernel(deep, M_, estim_params_, oo_, options_, bayestopt_) iter=iter+1; end if iter==1000 @@ -214,7 +215,7 @@ for b=fpar:B if run_smoother [dr,info,M_,oo_] =compute_decision_rules(M_,opts_local,oo_); - if ismember(info(1),[3,4]) + if ismember(info(1),[3,4]) opts_local.qz_criterium = 1 + (opts_local.qz_criterium-1)*10; %loosen tolerance, useful for big models where BK conditions were tested on restricted state space [dr,info,M_,oo_] =compute_decision_rules(M_,opts_local,oo_); end diff --git a/matlab/prior_sampler.m b/matlab/prior_sampler.m index 75b1896fd..1303cbcf0 100644 --- a/matlab/prior_sampler.m +++ b/matlab/prior_sampler.m @@ -13,7 +13,7 @@ function results = prior_sampler(drsave,M_,bayestopt_,options_,oo_,estim_params_ % SPECIAL REQUIREMENTS % none -% Copyright © 2009-2018 Dynare Team +% Copyright © 2009-2023 Dynare Team % % This file is part of Dynare. % @@ -31,7 +31,7 @@ function results = prior_sampler(drsave,M_,bayestopt_,options_,oo_,estim_params_ % along with Dynare. If not, see . % Initialization. -prior_draw(bayestopt_, options_.prior_trunc); +Prior = dprior(bayestopt_, options_.prior_trunc); PriorDirectoryName = CheckPath('prior/draws',M_.dname); work = ~drsave; iteration = 0; @@ -94,10 +94,10 @@ while iteration < NumberOfSimulations dyn_waitbar(iteration/NumberOfSimulations,hh,'Please wait. Prior sampler...'); end loop_indx = loop_indx+1; - params = prior_draw(); - M_ = set_all_parameters(params,estim_params_,M_); - [T,R,~,INFO,M_,oo_] = dynare_resolve(M_,options_,oo_,'restrict'); - dr=oo_.dr; + params = Prior.draw(); + M_ = set_all_parameters(params, estim_params_, M_); + [T, R, ~, INFO, M_, oo_] = dynare_resolve(M_, options_, oo_, 'restrict'); + dr = oo_.dr; if ~INFO(1) INFO=endogenous_prior_restrictions(T,R,M_,options_,oo_); end @@ -200,4 +200,4 @@ m1 = m0 + (newobs'-m0)/iter; qq = m1*m1'; s1 = s0 + ( (newobs'*newobs-qq-s0) + (iter-1)*(m0*m0'-qq') )/iter; mu = m1; -sigma = s1; \ No newline at end of file +sigma = s1; diff --git a/matlab/slice_sampler.m b/matlab/slice_sampler.m index a50c12ab6..73fc6c347 100644 --- a/matlab/slice_sampler.m +++ b/matlab/slice_sampler.m @@ -24,7 +24,7 @@ function [theta, fxsim, neval] = slice_sampler(objective_function,theta,thetapri % SPECIAL REQUIREMENTS % none -% Copyright © 2015-2017 Dynare Team +% Copyright © 2015-2023 Dynare Team % % This file is part of Dynare. % @@ -58,6 +58,8 @@ neval = zeros(npar,1); % % % fname0=fname; fname = [ int2str(sampler_options.curr_block)]; +Prior = dprior(varargin{6},varargin{3}.prior_trunc); + it=0; while it= thetaprior(:,1)) && all(theta(:) <= thetaprior(:,2)) + theta = Prior.draw(); + if all(theta >= thetaprior(:,1)) && all(theta <= thetaprior(:,2)) fxold = -feval(objective_function,theta,varargin{:}); end end