From b4941c02d323fe1321cc56f90d80d52c50df5220 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Fri, 24 Apr 2015 18:12:46 +0200 Subject: [PATCH] Cosmetic Fixes to Metropolis-Hastings routines - Adds comments and headers and fixes typos in previous ones - Make naming in random_walk_metropolis_hastings_core.m more intuitive --- matlab/metropolis_hastings_initialization.m | 41 +++-- matlab/prior_posterior_statistics.m | 8 +- matlab/random_walk_metropolis_hastings.m | 45 ++--- matlab/random_walk_metropolis_hastings_core.m | 165 ++++++++---------- 4 files changed, 130 insertions(+), 129 deletions(-) diff --git a/matlab/metropolis_hastings_initialization.m b/matlab/metropolis_hastings_initialization.m index e8f4d0ef1..4650830e5 100644 --- a/matlab/metropolis_hastings_initialization.m +++ b/matlab/metropolis_hastings_initialization.m @@ -1,7 +1,7 @@ function [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ... metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_) -%function [ ix2, ilogpo2, ModelName, MhDirectoryName, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = -% metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds, dataset_,dataset_info,,options_,M_,estim_params_,bayestopt_,oo_) +%function [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ... +% metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_) % Metropolis-Hastings initialization. % % INPUTS @@ -11,6 +11,7 @@ function [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, % o vv [double] (p*p) matrix, posterior covariance matrix (at the mode). % o mh_bounds [double] (p*2) matrix defining lower and upper bounds for the parameters. % o dataset_ data structure +% o dataset_info dataset info structure % o options_ options structure % o M_ model structure % o estim_params_ estimated parameters structure @@ -18,12 +19,25 @@ function [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, % o oo_ outputs structure % % OUTPUTS -% None +% o ix2 [double] (nblck*npar) vector of starting points for different chains +% o ilogpo2 [double] (nblck*1) vector of initial posterior values for different chains +% o ModelName [string] name of the mod-file +% o MetropolisFolder [string] path to the Metropolis subfolder +% o fblck [scalar] number of the first MH chain to be run (not equal to 1 in case of recovery) +% o fline [double] (nblck*1) vector of first draw in each chain (not equal to 1 in case of recovery) +% o npar [scalar] number of parameters estimated +% o nblck [scalar] Number of MCM chains requested +% o nruns [double] (nblck*1) number of draws in each chain +% o NewFile [scalar] (nblck*1) vector storing the number +% of the first MH-file to created for each chain when saving additional +% draws +% o MAX_nruns [scalar] maximum number of draws in each MH-file on the harddisk +% o d [double] (p*p) matrix, Cholesky decomposition of the posterior covariance matrix (at the mode). % % SPECIAL REQUIREMENTS % None. -% Copyright (C) 2006-2013 Dynare Team +% Copyright (C) 2006-2015 Dynare Team % % This file is part of Dynare. % @@ -40,6 +54,7 @@ function [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . +%Initialize outputs ix2 = []; ilogpo2 = []; ModelName = []; @@ -68,7 +83,7 @@ MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8); d = chol(vv); if ~options_.load_mh_file && ~options_.mh_recover - % Here we start a new metropolis-hastings, previous draws are discarded. + % Here we start a new Metropolis-Hastings, previous draws are discarded. if nblck > 1 disp('Estimation::mcmc: Multiple chains mode.') else @@ -80,7 +95,7 @@ if ~options_.load_mh_file && ~options_.mh_recover delete([BaseName '_mh*_blck*.mat']); disp('Estimation::mcmc: Old mh-files successfully erased!') end - % Delete old metropolis log file. + % Delete old Metropolis log file. file = dir([ MetropolisFolder '/metropolis.log']); if length(file) delete([ MetropolisFolder '/metropolis.log']); @@ -128,11 +143,11 @@ if ~options_.load_mh_file && ~options_.mh_recover if init_iter > 100 && validate == 0 disp(['Estimation::mcmc: I couldn''t get a valid initial value in 100 trials.']) if options_.nointeractive - disp(['Estimation::mcmc: I reduce mh_init_scale by ten percent:']) + disp(['Estimation::mcmc: I reduce mh_init_scale by 10 percent:']) options_.mh_init_scale = .9*options_.mh_init_scale; disp(sprintf('Estimation::mcmc: Parameter mh_init_scale is now equal to %f.',options_.mh_init_scale)) else - disp(['Estimation::mcmc: You should Reduce mh_init_scale...']) + disp(['Estimation::mcmc: You should reduce mh_init_scale...']) disp(sprintf('Estimation::mcmc: Parameter mh_init_scale is equal to %f.',options_.mh_init_scale)) options_.mh_init_scale = input('Estimation::mcmc: Enter a new value... '); end @@ -217,7 +232,7 @@ if ~options_.load_mh_file && ~options_.mh_recover fclose(fidlog); elseif options_.load_mh_file && ~options_.mh_recover % Here we consider previous mh files (previous mh did not crash). - disp('Estimation::mcmc: I am loading past metropolis-hastings simulations...') + disp('Estimation::mcmc: I am loading past Metropolis-Hastings simulations...') load_last_mh_history_file(MetropolisFolder, ModelName); mh_files = dir([ MetropolisFolder filesep ModelName '_mh*.mat']); if ~length(mh_files) @@ -238,7 +253,7 @@ elseif options_.load_mh_file && ~options_.mh_recover nblck = past_number_of_blocks; options_.mh_nblck = nblck; end - % I read the last line of the last mh-file for initialization of the new metropolis-hastings simulations: + % I read the last line of the last mh-file for initialization of the new Metropolis-Hastings simulations: LastFileNumber = record.LastFileNumber; LastLineNumber = record.LastLineNumber; if LastLineNumber < MAX_nruns @@ -252,7 +267,7 @@ elseif options_.load_mh_file && ~options_.mh_recover ix2 = record.LastParameters; fblck = 1; NumberOfPreviousSimulations = sum(record.MhDraws(:,1),1); - fprintf('Estimation::mcmc: I am writting a new mh-history file... '); + fprintf('Estimation::mcmc: I am writing a new mh-history file... '); record.MhDraws = [record.MhDraws;zeros(1,3)]; NumberOfDrawsWrittenInThePastLastFile = MAX_nruns - LastLineNumber; NumberOfDrawsToBeSaved = nruns(1) - NumberOfDrawsWrittenInThePastLastFile; @@ -318,7 +333,7 @@ elseif options_.mh_recover disp('Estimation::mcmc: It appears that you don''t need to use the mh_recover option!') disp(' You have to edit the mod file and remove the mh_recover option') disp(' in the estimation command') - error() + error('Estimation::mcmc: mh_recover option not required!') end % I count the number of saved mh files per block. NumberOfMhFilesPerBlock = zeros(nblck,1); @@ -339,7 +354,7 @@ elseif options_.mh_recover end % How many mh-files are saved in this block? NumberOfSavedMhFilesInTheCrashedBlck = NumberOfMhFilesPerBlock(fblck); - % Correct the number of saved mh files if the crashed metropolis was not the first session (so + % Correct the number of saved mh files if the crashed Metropolis was not the first session (so % that NumberOfSavedMhFilesInTheCrashedBlck is the number of saved mh files in the crashed chain % of the current session). if OldMh diff --git a/matlab/prior_posterior_statistics.m b/matlab/prior_posterior_statistics.m index 4908d6b64..a2f43343e 100644 --- a/matlab/prior_posterior_statistics.m +++ b/matlab/prior_posterior_statistics.m @@ -14,11 +14,11 @@ function prior_posterior_statistics(type,dataset,dataset_info) % % SPECIAL REQUIREMENTS % none - +% % PARALLEL CONTEXT -% See the comments random_walk_metropolis_hastings.m funtion. - - +% See the comments in the random_walk_metropolis_hastings.m funtion. +% +% % Copyright (C) 2005-2013 Dynare Team % % This file is part of Dynare. diff --git a/matlab/random_walk_metropolis_hastings.m b/matlab/random_walk_metropolis_hastings.m index c2a49e7a2..2092922ac 100644 --- a/matlab/random_walk_metropolis_hastings.m +++ b/matlab/random_walk_metropolis_hastings.m @@ -1,14 +1,17 @@ function random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_) -%function record=random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_) -% Random walk Metropolis-Hastings algorithm. +% function random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_) +% Random Walk Metropolis-Hastings algorithm. % % INPUTS % o TargetFun [char] string specifying the name of the objective % function (posterior kernel). +% o ProposalFun [char] string specifying the name of the proposal +% density % o xparam1 [double] (p*1) vector of parameters to be estimated (initial values). % o vv [double] (p*p) matrix, posterior covariance matrix (at the mode). % o mh_bounds [double] (p*2) matrix defining lower and upper bounds for the parameters. % o dataset_ data structure +% o dataset_info dataset info structure % o options_ options structure % o M_ model structure % o estim_params_ estimated parameters structure @@ -16,27 +19,27 @@ function random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bou % o oo_ outputs structure % % ALGORITHM -% Metropolis-Hastings. +% Random-Walk Metropolis-Hastings. % % SPECIAL REQUIREMENTS % None. % % PARALLEL CONTEXT % The most computationally intensive part of this function may be executed -% in parallel. The code sutable to be executed in -% parallel on multi core or cluster machine (in general a 'for' cycle), -% is removed from this function and placed in random_walk_metropolis_hastings_core.m funtion. -% Then the DYNARE parallel package contain a set of pairs matlab functions that can be executed in -% parallel and called name_function.m and name_function_core.m. -% In addition in parallel package we have second set of functions used -% to manage the parallel computation. +% in parallel. The code suitable to be executed in +% parallel on multi core or cluster machine (in general a 'for' cycle) +% has been removed from this function and been placed in the random_walk_metropolis_hastings_core.m funtion. +% +% The DYNARE parallel packages comprise a i) set of pairs of Matlab functions that can be executed in +% parallel and called name_function.m and name_function_core.m and ii) a second set of functions used +% to manage the parallel computations. % -% This function was the first function to be parallelized, later other +% This function was the first function to be parallelized. Later, other % functions have been parallelized using the same methodology. % Then the comments write here can be used for all the other pairs of -% parallel functions and also for management funtions. +% parallel functions and also for management functions. -% Copyright (C) 2006-2013 Dynare Team +% Copyright (C) 2006-2015 Dynare Team % % This file is part of Dynare. % @@ -54,7 +57,7 @@ function random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bou % along with Dynare. If not, see . -% In Metropolis, we set penalty to Inf to as to reject all parameter sets triggering error in target density computation +% In Metropolis, we set penalty to Inf so as to reject all parameter sets triggering an error during target density computation global objective_function_penalty_base objective_function_penalty_base = Inf; @@ -73,16 +76,16 @@ load_last_mh_history_file(MetropolisFolder, ModelName); % First run in serial mode, and then comment the follow line. % save('recordSerial.mat','-struct', 'record'); -% For parllel runs after serial runs with the abobe line active. +% For parallel runs after serial runs with the abobe line active. % TempRecord=load('recordSerial.mat'); % record.Seeds=TempRecord.Seeds; % Snapshot of the current state of computing. It necessary for the parallel -% execution (i.e. to execute in a corretct way portion of code remotely or -% on many core). The mandatory variables for local/remote parallel -% computing are stored in localVars struct. +% execution (i.e. to execute in a corretct way a portion of code remotely or +% on many cores). The mandatory variables for local/remote parallel +% computing are stored in the localVars struct. localVars = struct('TargetFun', TargetFun, ... 'ProposalFun', ProposalFun, ... @@ -110,9 +113,8 @@ localVars = struct('TargetFun', TargetFun, ... 'varargin',[]); -% The user don't want to use parallel computing, or want to compute a -% single chain. In this cases Random walk Metropolis-Hastings algorithm is -% computed sequentially. +% User doesn't want to use parallel computing, or wants to compute a +% single chain compute Random walk Metropolis-Hastings algorithm sequentially. if isnumeric(options_.parallel) || (nblck-fblck)==0, fout = random_walk_metropolis_hastings_core(localVars, fblck, nblck, 0); @@ -152,6 +154,7 @@ NewFile = fout(1).NewFile; update_last_mh_history_file(MetropolisFolder, ModelName, record); +% Provide diagnostic output skipline() disp(['Estimation::mcmc: Number of mh files: ' int2str(NewFile(1)) ' per block.']) disp(['Estimation::mcmc: Total number of generated files: ' int2str(NewFile(1)*nblck) '.']) diff --git a/matlab/random_walk_metropolis_hastings_core.m b/matlab/random_walk_metropolis_hastings_core.m index 90c89e871..2c26de0e5 100644 --- a/matlab/random_walk_metropolis_hastings_core.m +++ b/matlab/random_walk_metropolis_hastings_core.m @@ -1,25 +1,25 @@ function myoutput = random_walk_metropolis_hastings_core(myinputs,fblck,nblck,whoiam, ThisMatlab) -% PARALLEL CONTEXT -% This function contain the most computationally intensive portion of code in -% random_walk_metropolis_hastings (the 'for xxx = fblck:nblck' loop). The branches in 'for' -% cycle and are completely independent than suitable to be executed in parallel way. +% function myoutput = random_walk_metropolis_hastings_core(myinputs,fblck,nblck,whoiam, ThisMatlab) +% Contains the most computationally intensive portion of code in +% random_walk_metropolis_hastings (the 'for xxx = fblck:nblck' loop). The branches in that 'for' +% cycle are completely independent to be suitable for parallel execution. % % INPUTS % o myimput [struc] The mandatory variables for local/remote % parallel computing obtained from random_walk_metropolis_hastings.m % function. % o fblck and nblck [integer] The Metropolis-Hastings chains. -% o whoiam [integer] In concurrent programming a modality to refer to the differents thread running in parallel is needed. +% o whoiam [integer] In concurrent programming a modality to refer to the different threads running in parallel is needed. % The integer whoaim is the integer that % allows us to distinguish between them. Then it is the index number of this CPU among all CPUs in the % cluster. % o ThisMatlab [integer] Allows us to distinguish between the -% 'main' matlab, the slave matlab worker, local matlab, remote matlab, +% 'main' Matlab, the slave Matlab worker, local Matlab, remote Matlab, % ... Then it is the index number of this slave machine in the cluster. % OUTPUTS % o myoutput [struc] -% If executed without parallel is the original output of 'for b = -% fblck:nblck' otherwise a portion of it computed on a specific core or +% If executed without parallel, this is the original output of 'for b = +% fblck:nblck'. Otherwise, it's a portion of it computed on a specific core or % remote machine. In this case: % record; % irun; @@ -31,23 +31,12 @@ function myoutput = random_walk_metropolis_hastings_core(myinputs,fblck,nblck,wh % % SPECIAL REQUIREMENTS. % None. - +% % PARALLEL CONTEXT -% The most computationally intensive part of this function may be executed -% in parallel. The code sutable to be executed in parallel on multi core or cluster machine, -% is removed from this function and placed in random_walk_metropolis_hastings_core.m funtion. -% Then the DYNARE parallel package contain a set of pairs matlab functios that can be executed in -% parallel and called name_function.m and name_function_core.m. -% In addition in the parallel package we have second set of functions used -% to manage the parallel computation. -% -% This function was the first function to be parallelized, later other -% functions have been parallelized using the same methodology. -% Then the comments write here can be used for all the other pairs of -% parallel functions and also for management funtions. +% See the comments in the random_walk_metropolis_hastings.m funtion. -% Copyright (C) 2006-2013 Dynare Team +% Copyright (C) 2006-2015 Dynare Team % % This file is part of Dynare. % @@ -74,11 +63,9 @@ end TargetFun=myinputs.TargetFun; ProposalFun=myinputs.ProposalFun; xparam1=myinputs.xparam1; -vv=myinputs.vv; mh_bounds=myinputs.mh_bounds; -ix2=myinputs.ix2; -ilogpo2=myinputs.ilogpo2; -ModelName=myinputs.ModelName; +last_draw=myinputs.ix2; +last_posterior=myinputs.ilogpo2; fline=myinputs.fline; npar=myinputs.npar; nruns=myinputs.nruns; @@ -94,11 +81,9 @@ estim_params_ = myinputs.estim_params_; options_ = myinputs.options_; M_ = myinputs.M_; oo_ = myinputs.oo_; -varargin=myinputs.varargin; % Necessary only for remote computing! if whoiam - Parallel=myinputs.Parallel; % initialize persistent variables in priordens() priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, bayestopt_.p3,bayestopt_.p4,1); end @@ -116,53 +101,51 @@ elseif strcmpi(ProposalFun,'rand_multivariate_student') end % -% NOW i run the (nblck-fblck+1) metropolis-hastings chains +% Now I run the (nblck-fblck+1) Metropolis-Hastings chains % proposal_covariance_Cholesky_decomposition = d*diag(bayestopt_.jscale); -jloop=0; +block_iter=0; -JSUM = 0; -for b = fblck:nblck, - jloop=jloop+1; +for curr_block = fblck:nblck, + block_iter=block_iter+1; try - % This will not work if the master uses a random generator not + % This will not work if the master uses a random number generator not % available in the slave (different Matlab version or - % Matlab/Octave cluster). Therefor the trap. + % Matlab/Octave cluster). Therefore the trap. % - % This set the random generator type (the seed is useless but - % needed by the function) + % Set the random number generator type (the seed is useless but needed by the function) set_dynare_seed(options_.DynareRandomStreams.algo, options_.DynareRandomStreams.seed); - % This set the state - set_dynare_random_generator_state(record.InitialSeeds(b).Unifor, record.InitialSeeds(b).Normal); + % Set the state of the RNG + set_dynare_random_generator_state(record.InitialSeeds(curr_block).Unifor, record.InitialSeeds(curr_block).Normal); catch - % If the state set by master is incompatible with the slave, we - % only reseed - set_dynare_seed(options_.DynareRandomStreams.seed+b); + % If the state set by master is incompatible with the slave, we only reseed + set_dynare_seed(options_.DynareRandomStreams.seed+curr_block); end - if (options_.load_mh_file~=0) && (fline(b)>1) && OpenOldFile(b) - load([BaseName '_mh' int2str(NewFile(b)) '_blck' int2str(b) '.mat']) - x2 = [x2;zeros(InitSizeArray(b)-fline(b)+1,npar)]; - logpo2 = [logpo2;zeros(InitSizeArray(b)-fline(b)+1,1)]; - OpenOldFile(b) = 0; + if (options_.load_mh_file~=0) && (fline(curr_block)>1) && OpenOldFile(curr_block) %load previous draws and likelihood + load([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat']) + x2 = [x2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,npar)]; + logpo2 = [logpo2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,1)]; + OpenOldFile(curr_block) = 0; else - x2 = zeros(InitSizeArray(b),npar); - logpo2 = zeros(InitSizeArray(b),1); + x2 = zeros(InitSizeArray(curr_block),npar); + logpo2 = zeros(InitSizeArray(curr_block),1); end + %Prepare waiting bars if whoiam - prc0=(b-fblck)/(nblck-fblck+1)*(isoctave || options_.console_mode); - hh = dyn_waitbar({prc0,whoiam,options_.parallel(ThisMatlab)},['MH (' int2str(b) '/' int2str(options_.mh_nblck) ')...']); + prc0=(curr_block-fblck)/(nblck-fblck+1)*(isoctave || options_.console_mode); + hh = dyn_waitbar({prc0,whoiam,options_.parallel(ThisMatlab)},['MH (' int2str(curr_block) '/' int2str(options_.mh_nblck) ')...']); else - hh = dyn_waitbar(0,['Metropolis-Hastings (' int2str(b) '/' int2str(options_.mh_nblck) ')...']); + hh = dyn_waitbar(0,['Metropolis-Hastings (' int2str(curr_block) '/' int2str(options_.mh_nblck) ')...']); set(hh,'Name','Metropolis-Hastings'); end - isux = 0; - jsux = 0; - irun = fline(b); - j = 1; - while j <= nruns(b) - par = feval(ProposalFun, ix2(b,:), proposal_covariance_Cholesky_decomposition, n); + accepted_draws_this_chain = 0; + accepted_draws_this_file = 0; + draw_index_current_file = fline(curr_block); %get location of first draw in current block + draw_iter = 1; + while draw_iter <= nruns(curr_block) + par = feval(ProposalFun, last_draw(curr_block,:), proposal_covariance_Cholesky_decomposition, n); if all( par(:) > mh_bounds.lb ) && all( par(:) < mh_bounds.ub ) try logpost = - feval(TargetFun, par(:),dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_); @@ -172,29 +155,29 @@ for b = fblck:nblck, else logpost = -inf; end - if (logpost > -inf) && (log(rand) < logpost-ilogpo2(b)) - x2(irun,:) = par; - ix2(b,:) = par; - logpo2(irun) = logpost; - ilogpo2(b) = logpost; - isux = isux + 1; - jsux = jsux + 1; + if (logpost > -inf) && (log(rand) < logpost-last_posterior(curr_block)) + x2(draw_index_current_file,:) = par; + last_draw(curr_block,:) = par; + logpo2(draw_index_current_file) = logpost; + last_posterior(curr_block) = logpost; + accepted_draws_this_chain = accepted_draws_this_chain + 1; + accepted_draws_this_file = accepted_draws_this_file + 1; else - x2(irun,:) = ix2(b,:); - logpo2(irun) = ilogpo2(b); + x2(draw_index_current_file,:) = last_draw(curr_block,:); + logpo2(draw_index_current_file) = last_posterior(curr_block); end - prtfrc = j/nruns(b); - if (mod(j, 3)==0 && ~whoiam) || (mod(j,50)==0 && whoiam) - dyn_waitbar(prtfrc,hh,[ 'MH (' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('Current acceptance ratio %4.3f', isux/j)]); + prtfrc = draw_iter/nruns(curr_block); + if (mod(draw_iter, 3)==0 && ~whoiam) || (mod(draw_iter,50)==0 && whoiam) + dyn_waitbar(prtfrc,hh,[ 'MH (' int2str(curr_block) '/' int2str(options_.mh_nblck) ') ' sprintf('Current acceptance ratio %4.3f', accepted_draws_this_chain/draw_iter)]); end - if (irun == InitSizeArray(b)) || (j == nruns(b)) % Now I save the simulations - save([BaseName '_mh' int2str(NewFile(b)) '_blck' int2str(b) '.mat'],'x2','logpo2'); + if (draw_index_current_file == InitSizeArray(curr_block)) || (draw_iter == nruns(curr_block)) % Now I save the simulations, either because the current file is full or the chain is done + save([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat'],'x2','logpo2'); fidlog = fopen([MetropolisFolder '/metropolis.log'],'a'); fprintf(fidlog,['\n']); - fprintf(fidlog,['%% Mh' int2str(NewFile(b)) 'Blck' int2str(b) ' (' datestr(now,0) ')\n']); + fprintf(fidlog,['%% Mh' int2str(NewFile(curr_block)) 'Blck' int2str(curr_block) ' (' datestr(now,0) ')\n']); fprintf(fidlog,' \n'); fprintf(fidlog,[' Number of simulations.: ' int2str(length(logpo2)) '\n']); - fprintf(fidlog,[' Acceptance ratio......: ' num2str(jsux/length(logpo2)) '\n']); + fprintf(fidlog,[' Acceptance ratio......: ' num2str(accepted_draws_this_file/length(logpo2)) '\n']); fprintf(fidlog,[' Posterior mean........:\n']); for i=1:length(x2(1,:)) fprintf(fidlog,[' params:' int2str(i) ': ' num2str(mean(x2(:,i))) '\n']); @@ -212,32 +195,32 @@ for b = fblck:nblck, fprintf(fidlog,[' log2po:' num2str(max(logpo2)) '\n']); fprintf(fidlog,' \n'); fclose(fidlog); - jsux = 0; - if j == nruns(b) % I record the last draw... - record.LastParameters(b,:) = x2(end,:); - record.LastLogPost(b) = logpo2(end); + accepted_draws_this_file = 0; + if draw_iter == nruns(curr_block) % I record the last draw... + record.LastParameters(curr_block,:) = x2(end,:); + record.LastLogPost(curr_block) = logpo2(end); end - % size of next file in chain b - InitSizeArray(b) = min(nruns(b)-j,MAX_nruns); + % size of next file in chain curr_block + InitSizeArray(curr_block) = min(nruns(curr_block)-draw_iter,MAX_nruns); % initialization of next file if necessary - if InitSizeArray(b) - x2 = zeros(InitSizeArray(b),npar); - logpo2 = zeros(InitSizeArray(b),1); - NewFile(b) = NewFile(b) + 1; - irun = 0; + if InitSizeArray(curr_block) + x2 = zeros(InitSizeArray(curr_block),npar); + logpo2 = zeros(InitSizeArray(curr_block),1); + NewFile(curr_block) = NewFile(curr_block) + 1; + draw_index_current_file = 0; end end - j=j+1; - irun = irun + 1; + draw_iter=draw_iter+1; + draw_index_current_file = draw_index_current_file + 1; end% End of the simulations for one mh-block. - record.AcceptanceRatio(b) = isux/j; + record.AcceptanceRatio(curr_block) = accepted_draws_this_chain/draw_iter; dyn_waitbar_close(hh); - [record.LastSeeds(b).Unifor, record.LastSeeds(b).Normal] = get_dynare_random_generator_state(); - OutputFileName(jloop,:) = {[MetropolisFolder,filesep], [ModelName '_mh*_blck' int2str(b) '.mat']}; + [record.LastSeeds(curr_block).Unifor, record.LastSeeds(curr_block).Normal] = get_dynare_random_generator_state(); + OutputFileName(block_iter,:) = {[MetropolisFolder,filesep], [ModelName '_mh*_blck' int2str(curr_block) '.mat']}; end% End of the loop over the mh-blocks. myoutput.record = record; -myoutput.irun = irun; +myoutput.irun = draw_index_current_file; myoutput.NewFile = NewFile; myoutput.OutputFileName = OutputFileName; \ No newline at end of file