From 3cfcc98ed7c0a698a9e81654d463615c07421d14 Mon Sep 17 00:00:00 2001
From: Johannes Pfeifer <jpfeifer@gmx.de>
Date: Thu, 3 Oct 2024 20:44:43 +0200
Subject: [PATCH] GetAllPosteriorDraws.m: make routine more efficient by
enabling to output full parameters matrix instead of column vectors
Also factorizes loading of record file to obtain MCMC info by moving it into GetAllPosteriorDraws.m
---
.../mcmc_diagnostics.m | 7 +-
matlab/estimation/GetAllPosteriorDraws.m | 91 +++++++++++++------
.../GetPosteriorParametersStatistics.m | 33 +------
.../estimation/mh_autocorrelation_function.m | 15 +--
.../estimation/prior_posterior_statistics.m | 18 +---
matlab/estimation/trace_plot.m | 23 +----
6 files changed, 82 insertions(+), 105 deletions(-)
diff --git a/matlab/convergence_diagnostics/mcmc_diagnostics.m b/matlab/convergence_diagnostics/mcmc_diagnostics.m
index 6872074fb..ad216be37 100644
--- a/matlab/convergence_diagnostics/mcmc_diagnostics.m
+++ b/matlab/convergence_diagnostics/mcmc_diagnostics.m
@@ -59,16 +59,14 @@ if issue_an_error_message
end
% compute inefficiency factor
-FirstLine = record.KeepedDraws.FirstLine;
-TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
-FirstMhFile = record.KeepedDraws.FirstMhFile;
NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
param_name = {};
param_name_tex = {};
Ifac=NaN(nblck,npar);
+AllDraws = GetAllPosteriorDraws(options_, M_.dname, M_.fname, 'all');
for jj = 1:npar
if options_.TeX
[par_name_temp, par_name_tex_temp] = get_the_name(jj, options_.TeX, M_,estim_params_, options_.varobs);
@@ -79,8 +77,7 @@ for jj = 1:npar
par_name_temp = get_the_name(jj, options_.TeX, M_, estim_params_, options_.varobs);
param_name = vertcat(param_name, par_name_temp);
end
- Draws = GetAllPosteriorDraws(options_, M_.dname, M_.fname, jj, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck);
- Draws = reshape(Draws, [NumberOfDraws nblck]);
+ Draws = reshape(AllDraws(:,jj), [NumberOfDraws nblck]);
Nc = min(1000, NumberOfDraws/2);
for ll = 1:nblck
Ifac(ll,jj) = mcmc_ifac(Draws(:,ll), Nc);
diff --git a/matlab/estimation/GetAllPosteriorDraws.m b/matlab/estimation/GetAllPosteriorDraws.m
index fabf390af..6c1fd0211 100644
--- a/matlab/estimation/GetAllPosteriorDraws.m
+++ b/matlab/estimation/GetAllPosteriorDraws.m
@@ -1,26 +1,22 @@
-function draws = GetAllPosteriorDraws(options_, dname, fname, column, FirstMhFile, FirstLine, TotalNumberOfMhFile, NumberOfDraws, nblcks, blck)
-
+function draws = GetAllPosteriorDraws(options_, dname, fname, column, FirstLine , FirstMhFile, blck)
+% function draws = GetAllPosteriorDraws(options_, dname, fname, column, FirstLine , FirstMhFile, blck)
% Gets all posterior draws.
%
% INPUTS
-% - options_ [struct] Dynare's options.
-% - dname [char] name of directory with results.
-% - fname [char] name of mod file.
-% - column [integer] scalar, parameter index.
-% - FirstMhFile [integer] scalar, first MH file.
-% - FirstLine [integer] scalar, first line in first MH file.
-% - TotalNumberOfMhFile [integer] scalar, total number of MH file.
-% - NumberOfDraws [integer] scalar, number of posterior draws.
-% - nblcks [integer] scalar, total number of blocks.
-% - blck: [integer] scalar, desired block to read.
+% - options_ [struct] Dynare's options.
+% - dname [char] name of directory with results.
+% - fname [char] name of mod file.
+% - column [integer or string] scalar, parameter index
+% 'all' if all parameters are
+% requested
+% - FirstMhFile [integer] optional scalar, first MH file to be used
+% - FirstLine [integer] optional scalar, first line in first MH file to be used
+% - blck: [integer] scalar, desired block to read (optional)
%
% OUTPUTS
% - draws: [double] NumberOfDraws×1 vector, draws from posterior distribution.
-%
-% REMARKS
-% Only the first and third input arguments are required for SMC samplers.
-% Copyright © 2005-2023 Dynare Team
+% Copyright © 2005-2024 Dynare Team
%
% This file is part of Dynare.
%
@@ -37,23 +33,35 @@ function draws = GetAllPosteriorDraws(options_, dname, fname, column, FirstMhFil
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+if isnumeric(column) && column<0
+ error('GetAllPosteriorDraws:: column cannot be negative');
+end
+
if issmc(options_)
if ishssmc(options_)
% Load draws from the posterior distribution
pfiles = dir(sprintf('%s/hssmc/particles-*.mat', dname));
posterior = load(sprintf('%s/hssmc/particles-%u-%u.mat', dname, length(pfiles), length(pfiles)));
- if column==0
- draws = posterior.tlogpostkernel;
+ if column>0 || strcmp(column,'all')
+ if strcmp(column,'all')
+ draws = transpose(posterior.particles);
+ else
+ draws = transpose(posterior.particles(column,:));
+ end
else
- draws = transpose(posterior.particles(column,:));
+ draws = posterior.tlogpostkernel;
end
elseif isdime(options_)
posterior = load(sprintf('%s%s%s%schains.mat', dname, filesep(), 'dime', filesep()));
tune = posterior.tune;
chains = posterior.chains(end-tune:end,:,:);
- if column>0
+ if column>0 || strcmp(column,'all')
chains = reshape(chains, [], size(chains, 3));
- draws = chains(:,column);
+ if strcmp(column,'all')
+ draws = chains;
+ else
+ draws = chains(:,column);
+ end
else
draws = posterior.lprobs;
end
@@ -61,11 +69,26 @@ if issmc(options_)
error('GetAllPosteriorDraws:: case should not happen. Please contact the developers')
end
else
+ DirectoryName = CheckPath('metropolis',dname);
+ record=load_last_mh_history_file(DirectoryName,fname);
+ if nargin<5 || isempty(FirstLine)
+ FirstLine = record.KeepedDraws.FirstLine;
+ end
+ if nargin<6 || isempty(FirstMhFile)
+ FirstMhFile = record.KeepedDraws.FirstMhFile;
+ end
+ TotalNumberOfMhFile = sum(record.MhDraws(:,2));
+ TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
+ NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
+ [nblcks, npar] = size(record.LastParameters);
iline = FirstLine;
linee = 1;
- DirectoryName = CheckPath('metropolis',dname);
- if nblcks>1 && nargin<10
- draws = zeros(NumberOfDraws*nblcks,1);
+ if nblcks>1 && nargin<7
+ if strcmp(column,'all')
+ draws = zeros(NumberOfDraws*nblcks,npar);
+ else
+ draws = zeros(NumberOfDraws*nblcks,1);
+ end
iline0=iline;
if column>0
for blck = 1:nblcks
@@ -73,7 +96,11 @@ else
for file = FirstMhFile:TotalNumberOfMhFile
load([DirectoryName '/' fname '_mh' int2str(file) '_blck' int2str(blck)],'x2')
NumberOfLines = size(x2(iline:end,:),1);
- draws(linee:linee+NumberOfLines-1) = x2(iline:end,column);
+ if strcmp(column,'all')
+ draws(linee:linee+NumberOfLines-1,:) = x2(iline:end,:);
+ else
+ draws(linee:linee+NumberOfLines-1) = x2(iline:end,column);
+ end
linee = linee+NumberOfLines;
iline = 1;
end
@@ -90,15 +117,25 @@ else
end
end
end
- else
+ else %read desired block
if nblcks==1
blck=1;
end
+ if strcmp(column,'all')
+ draws = zeros(NumberOfDraws,npar);
+ else
+ draws = zeros(NumberOfDraws,1);
+ end
+
if column>0
for file = FirstMhFile:TotalNumberOfMhFile
load([DirectoryName '/' fname '_mh' int2str(file) '_blck' int2str(blck)],'x2')
NumberOfLines = size(x2(iline:end,:),1);
- draws(linee:linee+NumberOfLines-1) = x2(iline:end,column);
+ if strcmp(column,'all')
+ draws(linee:linee+NumberOfLines-1,:) = x2(iline:end,:);
+ else
+ draws(linee:linee+NumberOfLines-1) = x2(iline:end,column);
+ end
linee = linee+NumberOfLines;
iline = 1;
end
diff --git a/matlab/estimation/GetPosteriorParametersStatistics.m b/matlab/estimation/GetPosteriorParametersStatistics.m
index 299ccdaca..208a37779 100644
--- a/matlab/estimation/GetPosteriorParametersStatistics.m
+++ b/matlab/estimation/GetPosteriorParametersStatistics.m
@@ -44,17 +44,7 @@ np = estim_params_.np ;
latexFolder = CheckPath('latex',M_.dname);
FileName = M_.fname;
-if ~issmc(options_)
- MetropolisFolder = CheckPath('metropolis',M_.dname);
- record=load_last_mh_history_file(MetropolisFolder,FileName);
- FirstLine = record.KeepedDraws.FirstLine;
- TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
- TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
- FirstMhFile = record.KeepedDraws.FirstMhFile;
- NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
- mh_nblck = size(record.LastParameters,1);
- clear record;
-end
+[~, ~, num_draws]=set_number_of_subdraws(M_,options_); %get number of draws
header_width = row_header_width(M_, estim_params_, bayestopt_);
hpd_interval=[num2str(options_.mh_conf_sig*100), '% HPD interval'];
@@ -67,31 +57,18 @@ skipline()
if ishssmc(options_)
dprintf('Log data density is %f.', oo_.MarginalDensity.hssmc);
- % Set function handle for GetAllPosteriorDraws
- getalldraws = @(i) GetAllPosteriorDraws(options_, M_.dname, [], i);
+ hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
elseif isdime(options_)
- getalldraws = @(i) GetAllPosteriorDraws(options_, M_.dname, [], i);
+ hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
else
if ~isfield(oo_,'MarginalDensity') || (issmc(options_) && ~isfield(oo_.MarginalDensity,'ModifiedHarmonicMean'))
[~, oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_);
end
fprintf('Log data density (Modified Harmonic Mean) is %f.', oo_.MarginalDensity.ModifiedHarmonicMean);
- % Set function handle for GetAllPosteriordraws
- getalldraws = @(i) GetAllPosteriorDraws(options_, M_.dname, M_.fname, i, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, mh_nblck);
-end
-
-if ishssmc(options_)
- num_draws = options_.posterior_sampler_options.hssmc.particles;
- hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
-elseif isdime(options_)
- nchain = options_.posterior_sampler_options.current_options.nchain;
- tune = options_.posterior_sampler_options.current_options.tune;
- num_draws = nchain*tune;
- hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
-else
- num_draws=NumberOfDraws*mh_nblck;
hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
end
+% Set function handle for GetAllPosteriorDraws
+getalldraws = @(i) GetAllPosteriorDraws(options_, M_.dname, M_.fname, i);
if hpd_draws<2
fprintf('posterior_moments: There are not enough draws computes to compute HPD Intervals. Skipping their computation.\n')
diff --git a/matlab/estimation/mh_autocorrelation_function.m b/matlab/estimation/mh_autocorrelation_function.m
index 953b542ef..04ef6fb2c 100644
--- a/matlab/estimation/mh_autocorrelation_function.m
+++ b/matlab/estimation/mh_autocorrelation_function.m
@@ -45,21 +45,8 @@ if isempty(column)
return
end
-% Get informations about the posterior draws:
-MetropolisFolder = CheckPath('metropolis',M_.dname);
-record=load_last_mh_history_file(MetropolisFolder, M_.fname);
-
-FirstMhFile = record.KeepedDraws.FirstMhFile;
-FirstLine = record.KeepedDraws.FirstLine; ifil = FirstLine;
-TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
-TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
-NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
-nblck = size(record.LastParameters,1);
-
-clear record;
-
% Get all the posterior draws:
-PosteriorDraws = GetAllPosteriorDraws(options_, M_.dname, M_.fname, column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck, blck);
+PosteriorDraws = GetAllPosteriorDraws(options_, M_.dname, M_.fname, column, [],[], blck);
% Compute the autocorrelation function:
[~,autocor] = sample_autocovariance(PosteriorDraws,options_.mh_autocorrelation_function_size);
diff --git a/matlab/estimation/prior_posterior_statistics.m b/matlab/estimation/prior_posterior_statistics.m
index 519a58296..0768a4d30 100644
--- a/matlab/estimation/prior_posterior_statistics.m
+++ b/matlab/estimation/prior_posterior_statistics.m
@@ -203,20 +203,12 @@ localVars.bayestopt_=bayestopt_;
if strcmpi(type,'posterior')
- record=load_last_mh_history_file(DirectoryName, M_.fname);
- [nblck, npar] = size(record.LastParameters);
- FirstMhFile = record.KeepedDraws.FirstMhFile;
- FirstLine = record.KeepedDraws.FirstLine;
- TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
- TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
- NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
- mh_nblck = options_.mh_nblck;
- if B==NumberOfDraws*mh_nblck
+ [~, ~, NumberOfDraws]=set_number_of_subdraws(M_,options_);
+
+ if B==NumberOfDraws
% we load all retained MH runs !
- logpost=GetAllPosteriorDraws(options_, M_.dname, M_.fname, 0, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck);
- for column=1:npar
- x(:,column) = GetAllPosteriorDraws(options_, M_.dname, M_.fname, column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws, nblck);
- end
+ logpost=GetAllPosteriorDraws(options_, M_.dname, M_.fname, 0);
+ x = GetAllPosteriorDraws(options_, M_.dname, M_.fname, 'all');
else
[x, logpost]=get_posterior_subsample(M_,options_,B);
end
diff --git a/matlab/estimation/trace_plot.m b/matlab/estimation/trace_plot.m
index a3bb401ab..2d0dd7c2d 100644
--- a/matlab/estimation/trace_plot.m
+++ b/matlab/estimation/trace_plot.m
@@ -53,17 +53,7 @@ if isempty(column)
end
if ~issmc(options_)
- % Get informations about the posterior draws:
- MetropolisFolder = CheckPath('metropolis',M_.dname);
- record=load_last_mh_history_file(MetropolisFolder, M_.fname);
-
- FirstMhFile = 1;
- FirstLine = 1;
- TotalNumberOfMhFiles = sum(record.MhDraws(:,2));
- TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
- [mh_nblck] = size(record.LastParameters,2);
- clear record;
-
+ options_.mh_drop=0; % locally, do not drop as we want all draws
n_nblocks_to_plot=length(blck);
else
if ishssmc(options_)
@@ -73,14 +63,11 @@ else
end
end
+[~, ~, TotalNumberOfMhDraws]=set_number_of_subdraws(M_,options_);
+
if n_nblocks_to_plot==1
% Get all the posterior draws:
- if ishssmc(options_)
- PosteriorDraws = GetAllPosteriorDraws(options_, M_.dname,[],column);
- TotalNumberOfMhDraws=length(PosteriorDraws);
- else
- PosteriorDraws = GetAllPosteriorDraws(options_, M_.dname,M_.fname,column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, TotalNumberOfMhDraws, mh_nblck, blck);
- end
+ PosteriorDraws = GetAllPosteriorDraws(options_, M_.dname,M_.fname,column, 1, 1, blck);
else
PosteriorDraws=NaN(TotalNumberOfMhDraws,n_nblocks_to_plot);
save_string='';
@@ -88,7 +75,7 @@ else
title_string_tex='';
end
for block_iter=1:n_nblocks_to_plot
- PosteriorDraws(:,block_iter) = GetAllPosteriorDraws(options_, M_.dname, M_.fname, column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, TotalNumberOfMhDraws, mh_nblck, blck(block_iter));
+ PosteriorDraws(:,block_iter) = GetAllPosteriorDraws(options_, M_.dname, M_.fname, column, 1, 1, blck(block_iter));
save_string=[save_string,'_',num2str(blck(block_iter))];
if options_.TeX
title_string_tex=[title_string_tex, ', ' num2str(blck(block_iter))];
--
GitLab