From ca5c714e298e114ebddeca2643af3b1ea80de8ba Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Sat, 4 Apr 2015 19:36:57 +0200 Subject: [PATCH] Replace global variables in CutSample.m and metropolis_draw.m by function arguments Makes them usable outside of their current environment --- matlab/CutSample.m | 16 +++++----- matlab/dynare_estimation_1.m | 4 +-- matlab/metropolis_draw.m | 61 ++++++++++++++++++++++-------------- 3 files changed, 49 insertions(+), 32 deletions(-) diff --git a/matlab/CutSample.m b/matlab/CutSample.m index a3efbb01a..5ec0881a8 100644 --- a/matlab/CutSample.m +++ b/matlab/CutSample.m @@ -1,12 +1,14 @@ function CutSample(M_, options_, estim_params_) - -% function CutSample() -% Takes a subset from metropolis +% function CutSample(M_, options_, estim_params_) +% Takes a subset from metropolis draws by storing the required information +% like the first MH-file to be loaded and the first line in that file to be +% loaded into the record structure saved on harddisk into the +% _mh_history-file % % INPUTS -% options_ [structure] -% estim_params_ [structure] -% M_ [structure] +% M_ [structure] Dynare model structure +% options_ [structure] Dynare options structure +% estim_params_ [structure] Parameter structure % % OUTPUTS % none @@ -14,7 +16,7 @@ function CutSample(M_, options_, estim_params_) % SPECIAL REQUIREMENTS % none -% Copyright (C) 2005-2012 Dynare Team +% Copyright (C) 2005-2015 Dynare Team % % This file is part of Dynare. % diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index a813fbc0a..58282eba5 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -450,7 +450,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ... if ~options_.nodiagnostic && options_.mh_replic>0 oo_= McMCDiagnostics(options_, estim_params_, M_,oo_); end - %% Here i discard first half of the draws: + %% Here I discard first mh_drop percent of the draws: CutSample(M_, options_, estim_params_); %% Estimation of the marginal density from the Mh draws: if options_.mh_replic @@ -467,7 +467,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ... else load([M_.fname '_results'],'oo_'); end - error_flag = metropolis_draw(1); + [error_flag,junk,options_]= metropolis_draw(1,options_,estim_params_,M_); if options_.bayesian_irf if error_flag error('Estimation::mcmc: I cannot compute the posterior IRFs!') diff --git a/matlab/metropolis_draw.m b/matlab/metropolis_draw.m index 771ba0a43..f0fcf49bd 100644 --- a/matlab/metropolis_draw.m +++ b/matlab/metropolis_draw.m @@ -1,18 +1,30 @@ -function [xparams, logpost]=metropolis_draw(init) +function [xparams, logpost, options_]=metropolis_draw(init,options_,estim_params_,M_) % function [xparams, logpost]=metropolis_draw(init) % Builds draws from metropolis % % INPUTS: -% init: scalar equal to 1 (first call) or 0 +% init: scalar equal to +% 1: first call to store the required information +% on files, lines, and chains to be read +% in persistent variables to make them available +% for future calls +% 0: load a parameter draw +% Additional Inputs required for initialization +% options_ [structure] Matlab's structure describing the options (initialized by dynare, see @ref{options_}). +% M_ [structure] Matlab's structure describing the Model (initialized by dynare, see @ref{M_}). +% estim_params_ [structure] Matlab's structure describing the estimated_parameters (initialized by dynare, see @ref{estim_params_}). % % OUTPUTS: -% xparams: vector of estimated parameters +% xparams: if init==0: vector of estimated parameters +% if init==1: error flaog % logpost: log of posterior density +% options_: [structure] Matlab's structure describing the options (initialized by dynare, see @ref{options_}). % % SPECIAL REQUIREMENTS -% none +% +% Requires CutSample to be run before in order to set up mh_history-file -% Copyright (C) 2003-2011 Dynare Team +% Copyright (C) 2003-2015 Dynare Team % % This file is part of Dynare. % @@ -29,22 +41,24 @@ function [xparams, logpost]=metropolis_draw(init) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -global options_ estim_params_ M_ persistent mh_nblck NumberOfDraws BaseName FirstLine FirstMhFile MAX_nruns xparams = 0; logpost = 0; if init + %get number of parameters nvx = estim_params_.nvx; nvn = estim_params_.nvn; ncx = estim_params_.ncx; ncn = estim_params_.ncn; np = estim_params_.np ; npar = nvx+nvn+ncx+ncn+np; + %get path of metropolis files MetropolisFolder = CheckPath('metropolis',M_.dname); FileName = M_.fname; BaseName = [MetropolisFolder filesep FileName]; + %load mh_history-file with info on what to load load_last_mh_history_file(MetropolisFolder, FileName); FirstMhFile = record.KeepedDraws.FirstMhFile; FirstLine = record.KeepedDraws.FirstLine; @@ -52,7 +66,7 @@ if init LastMhFile = TotalNumberOfMhFiles; TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws); - MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8); + MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8); %number of parameters plus posterior plus ? mh_nblck = options_.mh_nblck; % set sub_draws option if empty if isempty(options_.sub_draws) @@ -67,20 +81,21 @@ if init end end return -end +else %not initialization, return one draw + %get random draw from random chain + ChainNumber = ceil(rand*mh_nblck); + DrawNumber = ceil(rand*NumberOfDraws); -ChainNumber = ceil(rand*mh_nblck); -DrawNumber = ceil(rand*NumberOfDraws); - -if DrawNumber <= MAX_nruns-FirstLine+1 - MhFilNumber = FirstMhFile; - MhLine = FirstLine+DrawNumber-1; -else - DrawNumber = DrawNumber-(MAX_nruns-FirstLine+1); - MhFilNumber = FirstMhFile+ceil(DrawNumber/MAX_nruns); - MhLine = DrawNumber-(MhFilNumber-FirstMhFile-1)*MAX_nruns; -end - -load( [ BaseName '_mh' int2str(MhFilNumber) '_blck' int2str(ChainNumber) '.mat' ],'x2','logpo2'); -xparams = x2(MhLine,:); -logpost= logpo2(MhLine); \ No newline at end of file + if DrawNumber <= MAX_nruns-FirstLine+1 %draw in first file, needs to account for first line + MhFilNumber = FirstMhFile; + MhLine = FirstLine+DrawNumber-1; + else %draw in other file + DrawNumber = DrawNumber-(MAX_nruns-FirstLine+1); + MhFilNumber = FirstMhFile+ceil(DrawNumber/MAX_nruns); + MhLine = DrawNumber-(MhFilNumber-FirstMhFile-1)*MAX_nruns; + end + %load parameters and posterior + load( [ BaseName '_mh' int2str(MhFilNumber) '_blck' int2str(ChainNumber) '.mat' ],'x2','logpo2'); + xparams = x2(MhLine,:); + logpost= logpo2(MhLine); +end \ No newline at end of file