stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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
time-shift
Johannes Pfeifer 2015-04-24 18:12:46 +02:00
parent 4ffebeee2e
commit b4941c02d3
4 changed files with 130 additions and 129 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

41
matlab/metropolis_hastings_initialization.m
View File

@ -1,7 +1,7 @@
function [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ... 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(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 ] = %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(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
% Metropolis-Hastings initialization. % Metropolis-Hastings initialization.
% %
% INPUTS % 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 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 mh_bounds [double] (p*2) matrix defining lower and upper bounds for the parameters.
% o dataset_ data structure % o dataset_ data structure
% o dataset_info dataset info structure
% o options_ options structure % o options_ options structure
% o M_ model structure % o M_ model structure
% o estim_params_ estimated parameters structure % o estim_params_ estimated parameters structure
@ -18,12 +19,25 @@ function [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck,
% o oo_ outputs structure % o oo_ outputs structure
% %
% OUTPUTS % 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 % SPECIAL REQUIREMENTS
% None. % None.
% Copyright (C) 2006-2013 Dynare Team % Copyright (C) 2006-2015 Dynare Team
% %
% This file is part of Dynare. % 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 % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>. % along with Dynare. If not, see <http://www.gnu.org/licenses/>.
%Initialize outputs
ix2 = []; ix2 = [];
ilogpo2 = []; ilogpo2 = [];
ModelName = []; ModelName = [];
@ -68,7 +83,7 @@ MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
d = chol(vv); d = chol(vv);
if ~options_.load_mh_file && ~options_.mh_recover 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 if nblck > 1
disp('Estimation::mcmc: Multiple chains mode.') disp('Estimation::mcmc: Multiple chains mode.')
else else
@ -80,7 +95,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
delete([BaseName '_mh*_blck*.mat']); delete([BaseName '_mh*_blck*.mat']);
disp('Estimation::mcmc: Old mh-files successfully erased!') disp('Estimation::mcmc: Old mh-files successfully erased!')
end end
% Delete old metropolis log file. % Delete old Metropolis log file.
file = dir([ MetropolisFolder '/metropolis.log']); file = dir([ MetropolisFolder '/metropolis.log']);
if length(file) if length(file)
delete([ MetropolisFolder '/metropolis.log']); delete([ MetropolisFolder '/metropolis.log']);
@ -128,11 +143,11 @@ if ~options_.load_mh_file && ~options_.mh_recover
if init_iter > 100 && validate == 0 if init_iter > 100 && validate == 0
disp(['Estimation::mcmc: I couldn''t get a valid initial value in 100 trials.']) disp(['Estimation::mcmc: I couldn''t get a valid initial value in 100 trials.'])
if options_.nointeractive 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; 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)) disp(sprintf('Estimation::mcmc: Parameter mh_init_scale is now equal to %f.',options_.mh_init_scale))
else 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)) 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... '); options_.mh_init_scale = input('Estimation::mcmc: Enter a new value... ');
end end
@ -217,7 +232,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
fclose(fidlog); fclose(fidlog);
elseif options_.load_mh_file && ~options_.mh_recover elseif options_.load_mh_file && ~options_.mh_recover
% Here we consider previous mh files (previous mh did not crash). % 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); load_last_mh_history_file(MetropolisFolder, ModelName);
mh_files = dir([ MetropolisFolder filesep ModelName '_mh*.mat']); mh_files = dir([ MetropolisFolder filesep ModelName '_mh*.mat']);
if ~length(mh_files) if ~length(mh_files)
@ -238,7 +253,7 @@ elseif options_.load_mh_file && ~options_.mh_recover
nblck = past_number_of_blocks; nblck = past_number_of_blocks;
options_.mh_nblck = nblck; options_.mh_nblck = nblck;
end 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; LastFileNumber = record.LastFileNumber;
LastLineNumber = record.LastLineNumber; LastLineNumber = record.LastLineNumber;
if LastLineNumber < MAX_nruns if LastLineNumber < MAX_nruns
@ -252,7 +267,7 @@ elseif options_.load_mh_file && ~options_.mh_recover
ix2 = record.LastParameters; ix2 = record.LastParameters;
fblck = 1; fblck = 1;
NumberOfPreviousSimulations = sum(record.MhDraws(:,1),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)]; record.MhDraws = [record.MhDraws;zeros(1,3)];
NumberOfDrawsWrittenInThePastLastFile = MAX_nruns - LastLineNumber; NumberOfDrawsWrittenInThePastLastFile = MAX_nruns - LastLineNumber;
NumberOfDrawsToBeSaved = nruns(1) - NumberOfDrawsWrittenInThePastLastFile; 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('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(' You have to edit the mod file and remove the mh_recover option')
disp(' in the estimation command') disp(' in the estimation command')
error() error('Estimation::mcmc: mh_recover option not required!')
end end
% I count the number of saved mh files per block. % I count the number of saved mh files per block.
NumberOfMhFilesPerBlock = zeros(nblck,1); NumberOfMhFilesPerBlock = zeros(nblck,1);
@ -339,7 +354,7 @@ elseif options_.mh_recover
end end
% How many mh-files are saved in this block? % How many mh-files are saved in this block?
NumberOfSavedMhFilesInTheCrashedBlck = NumberOfMhFilesPerBlock(fblck); 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 % that NumberOfSavedMhFilesInTheCrashedBlck is the number of saved mh files in the crashed chain
% of the current session). % of the current session).
if OldMh if OldMh

8
matlab/prior_posterior_statistics.m
View File

@ -14,11 +14,11 @@ function prior_posterior_statistics(type,dataset,dataset_info)
% %
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
%
% PARALLEL CONTEXT % 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 % Copyright (C) 2005-2013 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.

45
matlab/random_walk_metropolis_hastings.m
View File

@ -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 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_) % 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. % Random Walk Metropolis-Hastings algorithm.
% %
% INPUTS % INPUTS
% o TargetFun [char] string specifying the name of the objective % o TargetFun [char] string specifying the name of the objective
% function (posterior kernel). % 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 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 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 mh_bounds [double] (p*2) matrix defining lower and upper bounds for the parameters.
% o dataset_ data structure % o dataset_ data structure
% o dataset_info dataset info structure
% o options_ options structure % o options_ options structure
% o M_ model structure % o M_ model structure
% o estim_params_ estimated parameters 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 % o oo_ outputs structure
% %
% ALGORITHM % ALGORITHM
% Metropolis-Hastings. % Random-Walk Metropolis-Hastings.
% %
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% None. % None.
% %
% PARALLEL CONTEXT % PARALLEL CONTEXT
% The most computationally intensive part of this function may be executed % The most computationally intensive part of this function may be executed
% in parallel. The code sutable to be executed in % in parallel. The code suitable to be executed in
% parallel on multi core or cluster machine (in general a 'for' cycle), % 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. % has been removed from this function and been placed in the 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. % The DYNARE parallel packages comprise a i) set of pairs of Matlab functions that can be executed in
% In addition in parallel package we have second set of functions used % parallel and called name_function.m and name_function_core.m and ii) a second set of functions used
% to manage the parallel computation. % 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. % functions have been parallelized using the same methodology.
% Then the comments write here can be used for all the other pairs of % 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. % 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 <http://www.gnu.org/licenses/>. % along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% 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 global objective_function_penalty_base
objective_function_penalty_base = Inf; 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. % First run in serial mode, and then comment the follow line.
% save('recordSerial.mat','-struct', 'record'); % 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'); % TempRecord=load('recordSerial.mat');
% record.Seeds=TempRecord.Seeds; % record.Seeds=TempRecord.Seeds;
% Snapshot of the current state of computing. It necessary for the parallel % 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 % execution (i.e. to execute in a corretct way a portion of code remotely or
% on many core). The mandatory variables for local/remote parallel % on many cores). The mandatory variables for local/remote parallel
% computing are stored in localVars struct. % computing are stored in the localVars struct.
localVars = struct('TargetFun', TargetFun, ... localVars = struct('TargetFun', TargetFun, ...
'ProposalFun', ProposalFun, ... 'ProposalFun', ProposalFun, ...
@ -110,9 +113,8 @@ localVars = struct('TargetFun', TargetFun, ...
'varargin',[]); 'varargin',[]);
% The user don't want to use parallel computing, or want to compute a % User doesn't want to use parallel computing, or wants to compute a
% single chain. In this cases Random walk Metropolis-Hastings algorithm is % single chain compute Random walk Metropolis-Hastings algorithm sequentially.
% computed sequentially.
if isnumeric(options_.parallel) || (nblck-fblck)==0, if isnumeric(options_.parallel) || (nblck-fblck)==0,
fout = random_walk_metropolis_hastings_core(localVars, fblck, nblck, 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); update_last_mh_history_file(MetropolisFolder, ModelName, record);
% Provide diagnostic output
skipline() skipline()
disp(['Estimation::mcmc: Number of mh files: ' int2str(NewFile(1)) ' per block.']) disp(['Estimation::mcmc: Number of mh files: ' int2str(NewFile(1)) ' per block.'])
disp(['Estimation::mcmc: Total number of generated files: ' int2str(NewFile(1)*nblck) '.']) disp(['Estimation::mcmc: Total number of generated files: ' int2str(NewFile(1)*nblck) '.'])

165
matlab/random_walk_metropolis_hastings_core.m
View File

@ -1,25 +1,25 @@
function myoutput = random_walk_metropolis_hastings_core(myinputs,fblck,nblck,whoiam, ThisMatlab) function myoutput = random_walk_metropolis_hastings_core(myinputs,fblck,nblck,whoiam, ThisMatlab)
% PARALLEL CONTEXT % function myoutput = random_walk_metropolis_hastings_core(myinputs,fblck,nblck,whoiam, ThisMatlab)
% This function contain the most computationally intensive portion of code in % Contains the most computationally intensive portion of code in
% random_walk_metropolis_hastings (the 'for xxx = fblck:nblck' loop). The branches in 'for' % random_walk_metropolis_hastings (the 'for xxx = fblck:nblck' loop). The branches in that 'for'
% cycle and are completely independent than suitable to be executed in parallel way. % cycle are completely independent to be suitable for parallel execution.
% %
% INPUTS % INPUTS
% o myimput [struc] The mandatory variables for local/remote % o myimput [struc] The mandatory variables for local/remote
% parallel computing obtained from random_walk_metropolis_hastings.m % parallel computing obtained from random_walk_metropolis_hastings.m
% function. % function.
% o fblck and nblck [integer] The Metropolis-Hastings chains. % 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 % 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 % allows us to distinguish between them. Then it is the index number of this CPU among all CPUs in the
% cluster. % cluster.
% o ThisMatlab [integer] Allows us to distinguish between the % 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. % ... Then it is the index number of this slave machine in the cluster.
% OUTPUTS % OUTPUTS
% o myoutput [struc] % o myoutput [struc]
% If executed without parallel is the original output of 'for b = % If executed without parallel, this is the original output of 'for b =
% fblck:nblck' otherwise a portion of it computed on a specific core or % fblck:nblck'. Otherwise, it's a portion of it computed on a specific core or
% remote machine. In this case: % remote machine. In this case:
% record; % record;
% irun; % irun;
@ -31,23 +31,12 @@ function myoutput = random_walk_metropolis_hastings_core(myinputs,fblck,nblck,wh
% %
% SPECIAL REQUIREMENTS. % SPECIAL REQUIREMENTS.
% None. % None.
%
% PARALLEL CONTEXT % PARALLEL CONTEXT
% The most computationally intensive part of this function may be executed % See the comments in the random_walk_metropolis_hastings.m funtion.
% 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.
% Copyright (C) 2006-2013 Dynare Team % Copyright (C) 2006-2015 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -74,11 +63,9 @@ end
TargetFun=myinputs.TargetFun; TargetFun=myinputs.TargetFun;
ProposalFun=myinputs.ProposalFun; ProposalFun=myinputs.ProposalFun;
xparam1=myinputs.xparam1; xparam1=myinputs.xparam1;
vv=myinputs.vv;
mh_bounds=myinputs.mh_bounds; mh_bounds=myinputs.mh_bounds;
ix2=myinputs.ix2; last_draw=myinputs.ix2;
ilogpo2=myinputs.ilogpo2; last_posterior=myinputs.ilogpo2;
ModelName=myinputs.ModelName;
fline=myinputs.fline; fline=myinputs.fline;
npar=myinputs.npar; npar=myinputs.npar;
nruns=myinputs.nruns; nruns=myinputs.nruns;
@ -94,11 +81,9 @@ estim_params_ = myinputs.estim_params_;
options_ = myinputs.options_; options_ = myinputs.options_;
M_ = myinputs.M_; M_ = myinputs.M_;
oo_ = myinputs.oo_; oo_ = myinputs.oo_;
varargin=myinputs.varargin;
% Necessary only for remote computing! % Necessary only for remote computing!
if whoiam if whoiam
Parallel=myinputs.Parallel;
% initialize persistent variables in priordens() % initialize persistent variables in priordens()
priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, bayestopt_.p3,bayestopt_.p4,1); priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, bayestopt_.p3,bayestopt_.p4,1);
end end
@ -116,53 +101,51 @@ elseif strcmpi(ProposalFun,'rand_multivariate_student')
end 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); proposal_covariance_Cholesky_decomposition = d*diag(bayestopt_.jscale);
jloop=0; block_iter=0;
JSUM = 0; for curr_block = fblck:nblck,
for b = fblck:nblck, block_iter=block_iter+1;
jloop=jloop+1;
try 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 % 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 % Set the random number generator type (the seed is useless but needed by the function)
% needed by the function)
set_dynare_seed(options_.DynareRandomStreams.algo, options_.DynareRandomStreams.seed); set_dynare_seed(options_.DynareRandomStreams.algo, options_.DynareRandomStreams.seed);
% This set the state % Set the state of the RNG
set_dynare_random_generator_state(record.InitialSeeds(b).Unifor, record.InitialSeeds(b).Normal); set_dynare_random_generator_state(record.InitialSeeds(curr_block).Unifor, record.InitialSeeds(curr_block).Normal);
catch catch
% If the state set by master is incompatible with the slave, we % If the state set by master is incompatible with the slave, we only reseed
% only reseed set_dynare_seed(options_.DynareRandomStreams.seed+curr_block);
set_dynare_seed(options_.DynareRandomStreams.seed+b);
end end
if (options_.load_mh_file~=0) && (fline(b)>1) && OpenOldFile(b) if (options_.load_mh_file~=0) && (fline(curr_block)>1) && OpenOldFile(curr_block) %load previous draws and likelihood
load([BaseName '_mh' int2str(NewFile(b)) '_blck' int2str(b) '.mat']) load([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat'])
x2 = [x2;zeros(InitSizeArray(b)-fline(b)+1,npar)]; x2 = [x2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,npar)];
logpo2 = [logpo2;zeros(InitSizeArray(b)-fline(b)+1,1)]; logpo2 = [logpo2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,1)];
OpenOldFile(b) = 0; OpenOldFile(curr_block) = 0;
else else
x2 = zeros(InitSizeArray(b),npar); x2 = zeros(InitSizeArray(curr_block),npar);
logpo2 = zeros(InitSizeArray(b),1); logpo2 = zeros(InitSizeArray(curr_block),1);
end end
%Prepare waiting bars
if whoiam if whoiam
prc0=(b-fblck)/(nblck-fblck+1)*(isoctave || options_.console_mode); prc0=(curr_block-fblck)/(nblck-fblck+1)*(isoctave || options_.console_mode);
hh = dyn_waitbar({prc0,whoiam,options_.parallel(ThisMatlab)},['MH (' int2str(b) '/' int2str(options_.mh_nblck) ')...']); hh = dyn_waitbar({prc0,whoiam,options_.parallel(ThisMatlab)},['MH (' int2str(curr_block) '/' int2str(options_.mh_nblck) ')...']);
else 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'); set(hh,'Name','Metropolis-Hastings');
end end
isux = 0; accepted_draws_this_chain = 0;
jsux = 0; accepted_draws_this_file = 0;
irun = fline(b); draw_index_current_file = fline(curr_block); %get location of first draw in current block
j = 1; draw_iter = 1;
while j <= nruns(b) while draw_iter <= nruns(curr_block)
par = feval(ProposalFun, ix2(b,:), proposal_covariance_Cholesky_decomposition, n); par = feval(ProposalFun, last_draw(curr_block,:), proposal_covariance_Cholesky_decomposition, n);
if all( par(:) > mh_bounds.lb ) && all( par(:) < mh_bounds.ub ) if all( par(:) > mh_bounds.lb ) && all( par(:) < mh_bounds.ub )
try try
logpost = - feval(TargetFun, par(:),dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_); logpost = - feval(TargetFun, par(:),dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
@ -172,29 +155,29 @@ for b = fblck:nblck,
else else
logpost = -inf; logpost = -inf;
end end
if (logpost > -inf) && (log(rand) < logpost-ilogpo2(b)) if (logpost > -inf) && (log(rand) < logpost-last_posterior(curr_block))
x2(irun,:) = par; x2(draw_index_current_file,:) = par;
ix2(b,:) = par; last_draw(curr_block,:) = par;
logpo2(irun) = logpost; logpo2(draw_index_current_file) = logpost;
ilogpo2(b) = logpost; last_posterior(curr_block) = logpost;
isux = isux + 1; accepted_draws_this_chain = accepted_draws_this_chain + 1;
jsux = jsux + 1; accepted_draws_this_file = accepted_draws_this_file + 1;
else else
x2(irun,:) = ix2(b,:); x2(draw_index_current_file,:) = last_draw(curr_block,:);
logpo2(irun) = ilogpo2(b); logpo2(draw_index_current_file) = last_posterior(curr_block);
end end
prtfrc = j/nruns(b); prtfrc = draw_iter/nruns(curr_block);
if (mod(j, 3)==0 && ~whoiam) || (mod(j,50)==0 && whoiam) if (mod(draw_iter, 3)==0 && ~whoiam) || (mod(draw_iter,50)==0 && whoiam)
dyn_waitbar(prtfrc,hh,[ 'MH (' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('Current acceptance ratio %4.3f', isux/j)]); 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 end
if (irun == InitSizeArray(b)) || (j == nruns(b)) % Now I save the simulations 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(b)) '_blck' int2str(b) '.mat'],'x2','logpo2'); save([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat'],'x2','logpo2');
fidlog = fopen([MetropolisFolder '/metropolis.log'],'a'); fidlog = fopen([MetropolisFolder '/metropolis.log'],'a');
fprintf(fidlog,['\n']); 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,' \n');
fprintf(fidlog,[' Number of simulations.: ' int2str(length(logpo2)) '\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']); fprintf(fidlog,[' Posterior mean........:\n']);
for i=1:length(x2(1,:)) for i=1:length(x2(1,:))
fprintf(fidlog,[' params:' int2str(i) ': ' num2str(mean(x2(:,i))) '\n']); 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,[' log2po:' num2str(max(logpo2)) '\n']);
fprintf(fidlog,' \n'); fprintf(fidlog,' \n');
fclose(fidlog); fclose(fidlog);
jsux = 0; accepted_draws_this_file = 0;
if j == nruns(b) % I record the last draw... if draw_iter == nruns(curr_block) % I record the last draw...
record.LastParameters(b,:) = x2(end,:); record.LastParameters(curr_block,:) = x2(end,:);
record.LastLogPost(b) = logpo2(end); record.LastLogPost(curr_block) = logpo2(end);
end end
% size of next file in chain b % size of next file in chain curr_block
InitSizeArray(b) = min(nruns(b)-j,MAX_nruns); InitSizeArray(curr_block) = min(nruns(curr_block)-draw_iter,MAX_nruns);
% initialization of next file if necessary % initialization of next file if necessary
if InitSizeArray(b) if InitSizeArray(curr_block)
x2 = zeros(InitSizeArray(b),npar); x2 = zeros(InitSizeArray(curr_block),npar);
logpo2 = zeros(InitSizeArray(b),1); logpo2 = zeros(InitSizeArray(curr_block),1);
NewFile(b) = NewFile(b) + 1; NewFile(curr_block) = NewFile(curr_block) + 1;
irun = 0; draw_index_current_file = 0;
end end
end end
j=j+1; draw_iter=draw_iter+1;
irun = irun + 1; draw_index_current_file = draw_index_current_file + 1;
end% End of the simulations for one mh-block. 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); dyn_waitbar_close(hh);
[record.LastSeeds(b).Unifor, record.LastSeeds(b).Normal] = get_dynare_random_generator_state(); [record.LastSeeds(curr_block).Unifor, record.LastSeeds(curr_block).Normal] = get_dynare_random_generator_state();
OutputFileName(jloop,:) = {[MetropolisFolder,filesep], [ModelName '_mh*_blck' int2str(b) '.mat']}; OutputFileName(block_iter,:) = {[MetropolisFolder,filesep], [ModelName '_mh*_blck' int2str(curr_block) '.mat']};
end% End of the loop over the mh-blocks. end% End of the loop over the mh-blocks.
myoutput.record = record; myoutput.record = record;
myoutput.irun = irun; myoutput.irun = draw_index_current_file;
myoutput.NewFile = NewFile; myoutput.NewFile = NewFile;
myoutput.OutputFileName = OutputFileName; myoutput.OutputFileName = OutputFileName;