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

Added routines that factor out posterior samplers into the iteration routine posterior_sampler_iteration.m. 

It generalizes the random_walk_metropolis_hastings.m structure but adds two new items:
- check_posterior_sampler_options.m which is a wrapper for method specific options (stored in options_.posterior_options)
- posterior_sampler_iteration.m which is a wrapper of the individual iteration of any posterior sampler.

It already embeds, with parallelization implementation:
- random_walk metropolis hastings
- independent metropolis hastings
- slice sampler

It also contains a quick fix for TARB, although the latter should be embedded in posterior_sampler_iteration.m as well.

Any new posterior sampler can be simply added to posterior_sampler_iteration.m, with  wrapper to specific iterations in check_posterior_sampler_options.m
time-shift
Marco Ratto 2015-09-25 12:54:52 +02:00 committed by Johannes Pfeifer
parent b1b901fc61
commit c30b6f20a6
5 changed files with 1030 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

116
matlab/check_posterior_sampler_options.m Normal file
View File

@ -0,0 +1,116 @@
function [posterior_sampler_options, options_] = check_posterior_sampler_options(posterior_sampler_options, options_)
% function posterior_sampler_options = check_posterior_sampler_options(posterior_sampler_options, options_)
% initialization of posterior samplers
%
% INPUTS
% posterior_sampler_options: posterior sampler options
% options_: structure storing the options
% OUTPUTS
% posterior_sampler_options: checked posterior sampler options
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2015 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 <http://www.gnu.org/licenses/>.
posterior_sampler_options.posterior_sampling_method = options_.posterior_sampling_method;
posterior_sampler_options.proposal_distribution = options_.proposal_distribution;
bounds = posterior_sampler_options.bounds;
invhess = posterior_sampler_options.invhess;
% here are all samplers requiring a proposal distribution
if ~strcmp(posterior_sampler_options.posterior_sampling_method,'slice')
if ~options_.cova_compute
error('I Cannot start the MCMC because the Hessian of the posterior kernel at the mode was not computed.')
end
if options_.load_mh_file && options_.use_mh_covariance_matrix,
invhess = compute_mh_covariance_matrix;
posterior_sampler_options.invhess = invhess;
end
posterior_sampler_options.parallel_bar_refresh_rate=50;
posterior_sampler_options.serial_bar_refresh_rate=3;
posterior_sampler_options.parallel_bar_title='MH';
posterior_sampler_options.serial_bar_title='Metropolis-Hastings';
posterior_sampler_options.save_tmp_file=1;
end
% check specific options for slice sampler
if strcmp(options_.posterior_sampling_method,'slice')
posterior_sampler_options.parallel_bar_refresh_rate=1;
posterior_sampler_options.serial_bar_refresh_rate=1;
posterior_sampler_options.parallel_bar_title='SLICE';
posterior_sampler_options.serial_bar_title='SLICE';
posterior_sampler_options.save_tmp_file=1;
posterior_sampler_options = set_default_option(posterior_sampler_options,'rotated',0);
posterior_sampler_options = set_default_option(posterior_sampler_options,'slice_initialize_with_mode',0);
posterior_sampler_options = set_default_option(posterior_sampler_options,'use_slice_covariance_matrix',0);
posterior_sampler_options = set_default_option(posterior_sampler_options,'WR',[]);
if ~isfield(posterior_sampler_options,'mode'),
posterior_sampler_options.mode = [];
else % multimodal case
posterior_sampler_options.rotated = 1;
end
posterior_sampler_options = set_default_option(posterior_sampler_options,'mode_files',[]);
posterior_sampler_options = set_default_option(posterior_sampler_options,'W1',0.8*(bounds.ub-bounds.lb));
if options_.load_mh_file,
posterior_sampler_options.slice_initialize_with_mode = 0;
end
if posterior_sampler_options.rotated,
if isempty(posterior_sampler_options.mode_files) && ~isfield(posterior_sampler_options,'mode'), % rotated unimodal
if ~options_.cova_compute && ~(options_.load_mh_file && options_.use_mh_covariance_matrix)
skipline()
disp('I cannot start rotated slice sampler because')
disp('there is no previous MCMC to load ')
disp('or the Hessian at the mode is not computed.')
error('Rotated slice cannot start')
end
if isempty(invhess)
error('oops! This error should not occur, please contact developers.')
end
if options_.load_mh_file && posterior_sampler_options.use_slice_covariance_matrix,
invhess = compute_mh_covariance_matrix;
posterior_sampler_options.invhess = invhess;
end
[V1 D]=eig(invhess);
posterior_sampler_options.V1=V1;
posterior_sampler_options.WR=sqrt(diag(D))*3;
end
end
if ~isempty(posterior_sampler_options.mode_files), % multimodal case
modes = posterior_sampler_options.mode_files; % these can be also mean files from previous parallel slice chains
for j=1:length( modes ),
load(modes{j}, 'xparam1')
mode(j).m=xparam1;
end
posterior_sampler_options.mode = mode;
posterior_sampler_options.rotated = 1;
posterior_sampler_options.WR=[];
end
options_.mh_posterior_mode_estimation = 0;
end

196
matlab/posterior_sampler.m Normal file
View File

@ -0,0 +1,196 @@
function posterior_sampler(TargetFun,ProposalFun,xparam1,sampler_options,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
% function posterior_sampler(TargetFun,ProposalFun,xparam1,sampler_options,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 sampler_options structure
% .invhess [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
% o bayestopt_ estimation options structure
% o oo_ outputs structure
%
% ALGORITHM
% Random-Walk Metropolis-Hastings.
%
% SPECIAL REQUIREMENTS
% None.
%
% PARALLEL CONTEXT
% The most computationally intensive part of this function may be executed
% 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
% 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 functions.
% Copyright (C) 2006-2015 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 <http://www.gnu.org/licenses/>.
% 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;
vv = sampler_options.invhess;
% Initialization of the random walk metropolis-hastings chains.
[ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ...
posterior_sampler_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
InitSizeArray = min([repmat(MAX_nruns,nblck,1) fline+nruns-1],[],2);
% Load last mh history file
load_last_mh_history_file(MetropolisFolder, ModelName);
% Only for test parallel results!!!
% To check the equivalence between parallel and serial computation!
% First run in serial mode, and then comment the follow line.
% save('recordSerial.mat','-struct', 'record');
% 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 a portion of code remotely or
% on many cores). The mandatory variables for local/remote parallel
% computing are stored in the localVars struct.
if options_.TaRB.use_TaRB
options_.silent_optimizer=1; %locally set optimizer to silent mode
end
localVars = struct('TargetFun', TargetFun, ...
'ProposalFun', ProposalFun, ...
'xparam1', xparam1, ...
'vv', vv, ...
'sampler_options', sampler_options, ...
'mh_bounds', mh_bounds, ...
'ix2', ix2, ...
'ilogpo2', ilogpo2, ...
'ModelName', ModelName, ...
'fline', fline, ...
'npar', npar, ...
'nruns', nruns, ...
'NewFile', NewFile, ...
'MAX_nruns', MAX_nruns, ...
'd', d, ...
'InitSizeArray',InitSizeArray, ...
'record', record, ...
'dataset_', dataset_, ...
'dataset_info', dataset_info, ...
'options_', options_, ...
'M_',M_, ...
'bayestopt_', bayestopt_, ...
'estim_params_', estim_params_, ...
'oo_', oo_,...
'varargin',[]);
% 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,
if options_.TaRB.use_TaRB
fout = TaRB_metropolis_hastings_core(localVars, fblck, nblck, 0);
else
fout = posterior_sampler_core(localVars, fblck, nblck, 0);
end
record = fout.record;
% Parallel in Local or remote machine.
else
% Global variables for parallel routines.
globalVars = struct();
% which files have to be copied to run remotely
NamFileInput(1,:) = {'',[ModelName '_static.m']};
NamFileInput(2,:) = {'',[ModelName '_dynamic.m']};
if options_.steadystate_flag,
NamFileInput(length(NamFileInput)+1,:)={'',[ModelName '_steadystate.m']};
end
if (options_.load_mh_file~=0) && any(fline>1) ,
NamFileInput(length(NamFileInput)+1,:)={[M_.dname '/metropolis/'],[ModelName '_mh' int2str(NewFile(1)) '_blck*.mat']};
end
if exist([ModelName '_optimal_mh_scale_parameter.mat'])
NamFileInput(length(NamFileInput)+1,:)={'',[ModelName '_optimal_mh_scale_parameter.mat']};
end
% from where to get back results
% NamFileOutput(1,:) = {[M_.dname,'/metropolis/'],'*.*'};
if options_.TaRB.use_TaRB
[fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'TaRB_metropolis_hastings_core', localVars, globalVars, options_.parallel_info);
else
[fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'posterior_sampler_core', localVars, globalVars, options_.parallel_info);
end
for j=1:totCPU,
offset = sum(nBlockPerCPU(1:j-1))+fblck-1;
record.LastLogPost(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastLogPost(offset+1:sum(nBlockPerCPU(1:j)));
record.LastParameters(offset+1:sum(nBlockPerCPU(1:j)),:)=fout(j).record.LastParameters(offset+1:sum(nBlockPerCPU(1:j)),:);
record.AcceptanceRatio(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.AcceptanceRatio(offset+1:sum(nBlockPerCPU(1:j)));
record.FunctionEvalPerIteration(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.FunctionEvalPerIteration(offset+1:sum(nBlockPerCPU(1:j)));
record.LastSeeds(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastSeeds(offset+1:sum(nBlockPerCPU(1:j)));
end
end
irun = fout(1).irun;
NewFile = fout(1).NewFile;
record.MCMCConcludedSuccessfully = 1; %set indicator for successful run
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) '.'])
disp(['Estimation::mcmc: Total number of iterations: ' int2str((NewFile(1)-1)*MAX_nruns+irun-1) '.'])
disp(['Estimation::mcmc: Current acceptance ratio per chain: '])
for i=1:nblck
if i<10
disp([' Chain ' num2str(i) ': ' num2str(100*record.AcceptanceRatio(i)) '%'])
else
disp([' Chain ' num2str(i) ': ' num2str(100*record.AcceptanceRatio(i)) '%'])
end
end
if max(record.FunctionEvalPerIteration)>1
disp(['Estimation::mcmc: Current function evaluations per iteration: '])
for i=1:nblck
if i<10
disp([' Chain ' num2str(i) ': ' num2str(record.FunctionEvalPerIteration(i))])
else
disp([' Chain ' num2str(i) ': ' num2str(record.FunctionEvalPerIteration(i))])
end
end
end

241
matlab/posterior_sampler_core.m Normal file
View File

@ -0,0 +1,241 @@
function myoutput = posterior_sampler_core(myinputs,fblck,nblck,whoiam, ThisMatlab)
% function myoutput = posterior_sampler_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 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,
% ... Then it is the index number of this slave machine in the cluster.
% OUTPUTS
% o myoutput [struc]
% 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;
% NewFile;
% OutputFileName
%
% ALGORITHM
% Portion of Metropolis-Hastings.
%
% SPECIAL REQUIREMENTS.
% None.
%
% PARALLEL CONTEXT
% See the comments in the random_walk_metropolis_hastings.m funtion.
% Copyright (C) 2006-2015 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 <http://www.gnu.org/licenses/>.
if nargin<4,
whoiam=0;
end
% reshape 'myinputs' for local computation.
% In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by:
TargetFun=myinputs.TargetFun;
ProposalFun=myinputs.ProposalFun;
xparam1=myinputs.xparam1;
mh_bounds=myinputs.mh_bounds;
last_draw=myinputs.ix2;
last_posterior=myinputs.ilogpo2;
fline=myinputs.fline;
npar=myinputs.npar;
nruns=myinputs.nruns;
NewFile=myinputs.NewFile;
MAX_nruns=myinputs.MAX_nruns;
sampler_options=myinputs.sampler_options;
d=myinputs.d;
InitSizeArray=myinputs.InitSizeArray;
record=myinputs.record;
dataset_ = myinputs.dataset_;
dataset_info = myinputs.dataset_info;
bayestopt_ = myinputs.bayestopt_;
estim_params_ = myinputs.estim_params_;
options_ = myinputs.options_;
M_ = myinputs.M_;
oo_ = myinputs.oo_;
% Necessary only for remote computing!
if whoiam
% initialize persistent variables in priordens()
priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, bayestopt_.p3,bayestopt_.p4,1);
end
MetropolisFolder = CheckPath('metropolis',M_.dname);
ModelName = M_.fname;
BaseName = [MetropolisFolder filesep ModelName];
save_tmp_file = sampler_options.save_tmp_file;
options_.lik_algo = 1;
OpenOldFile = ones(nblck,1);
if strcmpi(ProposalFun,'rand_multivariate_normal')
sampler_options.n = npar;
sampler_options.ProposalDensity = 'multivariate_normal_pdf';
elseif strcmpi(ProposalFun,'rand_multivariate_student')
sampler_options.n = options_.student_degrees_of_freedom;
sampler_options.ProposalDensity = 'multivariate_student_pdf';
end
%
% Now I run the (nblck-fblck+1) MCMC chains
%
sampler_options.xparam1 = xparam1;
sampler_options.proposal_covariance_Cholesky_decomposition = d*diag(bayestopt_.jscale);
block_iter=0;
for curr_block = fblck:nblck,
block_iter=block_iter+1;
try
% This will not work if the master uses a random number generator not
% available in the slave (different Matlab version or
% Matlab/Octave cluster). Therefore the trap.
%
% Set the random number generator type (the seed is useless but needed by the function)
set_dynare_seed(options_.DynareRandomStreams.algo, options_.DynareRandomStreams.seed);
% 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+curr_block);
end
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(curr_block),npar);
logpo2 = zeros(InitSizeArray(curr_block),1);
end
%Prepare waiting bars
if whoiam
refresh_rate = sampler_options.parallel_bar_refresh_rate;
bar_title = sampler_options.parallel_bar_title;
prc0=(curr_block-fblck)/(nblck-fblck+1)*(isoctave || options_.console_mode);
hh = dyn_waitbar({prc0,whoiam,options_.parallel(ThisMatlab)},[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ')...']);
else
refresh_rate = sampler_options.serial_bar_refresh_rate;
bar_title = sampler_options.serial_bar_title;
hh = dyn_waitbar(0,[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ')...']);
set(hh,'Name',bar_title);
end
accepted_draws_this_chain = 0;
accepted_draws_this_file = 0;
feval_this_chain = 0;
feval_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, logpost, accepted, neval] = posterior_sampler_iteration(TargetFun, last_draw(curr_block,:), last_posterior(curr_block), sampler_options,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
x2(draw_index_current_file,:) = par;
last_draw(curr_block,:) = par;
logpo2(draw_index_current_file) = logpost;
last_posterior(curr_block) = logpost;
feval_this_chain = feval_this_chain + sum(neval);
feval_this_file = feval_this_file + sum(neval);
accepted_draws_this_chain = accepted_draws_this_chain + accepted;
accepted_draws_this_file = accepted_draws_this_file + accepted;
prtfrc = draw_iter/nruns(curr_block);
if mod(draw_iter, refresh_rate)==0
if accepted_draws_this_chain/draw_iter==1 && sum(neval)>1
dyn_waitbar(prtfrc,hh,[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ') ' sprintf('Function eval per draw %4.3f', feval_this_chain/draw_iter)]);
else
dyn_waitbar(prtfrc,hh,[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ') ' sprintf('Current acceptance ratio %4.3f', accepted_draws_this_chain/draw_iter)]);
end
if save_tmp_file
save([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat'],'x2','logpo2');
end
end
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
if save_tmp_file,
delete([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat']);
end
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(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(accepted_draws_this_file/length(logpo2)) '\n']);
fprintf(fidlog,[' Feval per iteration...: ' num2str(feval_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']);
end
fprintf(fidlog,[' log2po:' num2str(mean(logpo2)) '\n']);
fprintf(fidlog,[' Minimum value.........:\n']);
for i=1:length(x2(1,:))
fprintf(fidlog,[' params:' int2str(i) ': ' num2str(min(x2(:,i))) '\n']);
end
fprintf(fidlog,[' log2po:' num2str(min(logpo2)) '\n']);
fprintf(fidlog,[' Maximum value.........:\n']);
for i=1:length(x2(1,:))
fprintf(fidlog,[' params:' int2str(i) ': ' num2str(max(x2(:,i))) '\n']);
end
fprintf(fidlog,[' log2po:' num2str(max(logpo2)) '\n']);
fprintf(fidlog,' \n');
fclose(fidlog);
accepted_draws_this_file = 0;
feval_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 curr_block
InitSizeArray(curr_block) = min(nruns(curr_block)-draw_iter,MAX_nruns);
% initialization of next file if necessary
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
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(curr_block) = accepted_draws_this_chain/(draw_iter-1);
record.FunctionEvalPerIteration(curr_block) = feval_this_chain/(draw_iter-1);
dyn_waitbar_close(hh);
[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 = draw_index_current_file;
myoutput.NewFile = NewFile;
myoutput.OutputFileName = OutputFileName;

384
matlab/posterior_sampler_initialization.m Normal file
View File

@ -0,0 +1,384 @@
function [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ...
posterior_sampler_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
% o TargetFun [char] string specifying the name of the objective
% function (posterior kernel).
% 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
% o bayestopt_ estimation options structure
% o oo_ outputs structure
%
% OUTPUTS
% 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-2015 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 <http://www.gnu.org/licenses/>.
%Initialize outputs
ix2 = [];
ilogpo2 = [];
ModelName = [];
MetropolisFolder = [];
fblck = [];
fline = [];
npar = [];
nblck = [];
nruns = [];
NewFile = [];
MAX_nruns = [];
d = [];
ModelName = M_.fname;
if ~isempty(M_.bvar)
ModelName = [ModelName '_bvar'];
end
MetropolisFolder = CheckPath('metropolis',M_.dname);
BaseName = [MetropolisFolder filesep ModelName];
nblck = options_.mh_nblck;
nruns = ones(nblck,1)*options_.mh_replic;
npar = length(xparam1);
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.
if nblck > 1
disp('Estimation::mcmc: Multiple chains mode.')
else
disp('Estimation::mcmc: One Chain mode.')
end
% Delete old mh files if any...
files = dir([BaseName '_mh*_blck*.mat']);
if length(files)
delete([BaseName '_mh*_blck*.mat']);
disp('Estimation::mcmc: Old mh-files successfully erased!')
end
% Delete old Metropolis log file.
file = dir([ MetropolisFolder '/metropolis.log']);
if length(file)
delete([ MetropolisFolder '/metropolis.log']);
disp('Estimation::mcmc: Old metropolis.log file successfully erased!')
disp('Estimation::mcmc: Creation of a new metropolis.log file.')
end
fidlog = fopen([MetropolisFolder '/metropolis.log'],'w');
fprintf(fidlog,'%% MH log file (Dynare).\n');
fprintf(fidlog,['%% ' datestr(now,0) '.\n']);
fprintf(fidlog,' \n\n');
fprintf(fidlog,'%% Session 1.\n');
fprintf(fidlog,' \n');
fprintf(fidlog,[' Number of blocks...............: ' int2str(nblck) '\n']);
fprintf(fidlog,[' Number of simulations per block: ' int2str(nruns(1)) '\n']);
fprintf(fidlog,' \n');
% Find initial values for the nblck chains...
if isempty(d),
prior_draw(1);
end
if nblck > 1% Case 1: multiple chains
fprintf(fidlog,[' Initial values of the parameters:\n']);
disp('Estimation::mcmc: Searching for initial values...')
ix2 = zeros(nblck,npar);
ilogpo2 = zeros(nblck,1);
for j=1:nblck
validate = 0;
init_iter = 0;
trial = 1;
while validate == 0 && trial <= 10
if isempty(d),
candidate = prior_draw(0);
else
candidate = rand_multivariate_normal( transpose(xparam1), d * options_.mh_init_scale, npar);
end
if all(candidate(:) > mh_bounds.lb) && all(candidate(:) < mh_bounds.ub)
ix2(j,:) = candidate;
ilogpo2(j) = - feval(TargetFun,ix2(j,:)',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
if ~isfinite(ilogpo2(j)) % if returned log-density is
% Inf or Nan (penalized value)
validate = 0;
else
fprintf(fidlog,[' Blck ' int2str(j) ':\n']);
for i=1:length(ix2(1,:))
fprintf(fidlog,[' params:' int2str(i) ': ' num2str(ix2(j,i)) '\n']);
end
fprintf(fidlog,[' logpo2: ' num2str(ilogpo2(j)) '\n']);
j = j+1;
validate = 1;
end
end
init_iter = init_iter + 1;
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 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(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
trial = trial+1;
end
end
if trial > 10 && ~validate
disp(['Estimation::mcmc: I''m unable to find a starting value for block ' int2str(j)])
return
end
end
fprintf(fidlog,' \n');
disp('Estimation::mcmc: Initial values found!')
skipline()
else% Case 2: one chain (we start from the posterior mode)
fprintf(fidlog,[' Initial values of the parameters:\n']);
candidate = transpose(xparam1(:));%
if all(candidate(:) > mh_bounds.lb) && all(candidate(:) < mh_bounds.ub)
ix2 = candidate;
ilogpo2 = - feval(TargetFun,ix2',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
disp('Estimation::mcmc: Initialization at the posterior mode.')
skipline()
fprintf(fidlog,[' Blck ' int2str(1) 'params:\n']);
for i=1:length(ix2(1,:))
fprintf(fidlog,[' ' int2str(i) ':' num2str(ix2(1,i)) '\n']);
end
fprintf(fidlog,[' Blck ' int2str(1) 'logpo2:' num2str(ilogpo2) '\n']);
else
disp('Estimation::mcmc: Initialization failed...')
disp('Estimation::mcmc: The posterior mode lies outside the prior bounds.')
return
end
fprintf(fidlog,' \n');
end
fprintf(fidlog,' \n');
fblck = 1;
fline = ones(nblck,1);
NewFile = ones(nblck,1);
% Delete the mh-history files
delete_mh_history_files(MetropolisFolder, ModelName);
% Create a new record structure
fprintf(['Estimation::mcmc: Write details about the MCMC... ']);
AnticipatedNumberOfFiles = ceil(nruns(1)/MAX_nruns);
AnticipatedNumberOfLinesInTheLastFile = nruns(1) - (AnticipatedNumberOfFiles-1)*MAX_nruns;
record.Sampler = options_.posterior_sampling_method;
record.Nblck = nblck;
record.MhDraws = zeros(1,3);
record.MhDraws(1,1) = nruns(1);
record.MhDraws(1,2) = AnticipatedNumberOfFiles;
record.MhDraws(1,3) = AnticipatedNumberOfLinesInTheLastFile;
record.AcceptanceRatio = zeros(1,nblck);
record.FunctionEvalPerIteration = NaN(1,nblck);
for j=1:nblck
% 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);
% record.Seeds keeps a vector of the random generator state and not the scalar seed despite its name
[record.InitialSeeds(j).Unifor,record.InitialSeeds(j).Normal] = get_dynare_random_generator_state();
end
record.InitialParameters = ix2;
record.InitialLogPost = ilogpo2;
record.LastParameters = zeros(nblck,npar);
record.LastLogPost = zeros(nblck,1);
record.LastFileNumber = AnticipatedNumberOfFiles ;
record.LastLineNumber = AnticipatedNumberOfLinesInTheLastFile;
record.MCMCConcludedSuccessfully = 0;
fprintf('Ok!\n');
id = write_mh_history_file(MetropolisFolder, ModelName, record);
disp(['Estimation::mcmc: Details about the MCMC are available in ' BaseName '_mh_history_' num2str(id) '.mat'])
skipline()
fprintf(fidlog,[' CREATION OF THE MH HISTORY FILE!\n\n']);
fprintf(fidlog,[' Expected number of files per block.......: ' int2str(AnticipatedNumberOfFiles) '.\n']);
fprintf(fidlog,[' Expected number of lines in the last file: ' int2str(AnticipatedNumberOfLinesInTheLastFile) '.\n']);
fprintf(fidlog,['\n']);
for j = 1:nblck,
fprintf(fidlog,[' Initial state of the Gaussian random number generator for chain number ',int2str(j),':\n']);
for i=1:length(record.InitialSeeds(j).Normal)
fprintf(fidlog,[' ' num2str(record.InitialSeeds(j).Normal(i)') '\n']);
end
fprintf(fidlog,[' Initial state of the Uniform random number generator for chain number ',int2str(j),':\n']);
for i=1:length(record.InitialSeeds(j).Unifor)
fprintf(fidlog,[' ' num2str(record.InitialSeeds(j).Unifor(i)') '\n']);
end
end,
fprintf(fidlog,' \n');
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...')
load_last_mh_history_file(MetropolisFolder, ModelName);
if ~isnan(record.MCMCConcludedSuccessfully) && ~record.MCMCConcludedSuccessfully
error('Estimation::mcmc: You are trying to load an MCMC that did not finish successfully. Please use mh_recover.')
end
record.MCMCConcludedSuccessfully=0; %reset indicator for this run
mh_files = dir([ MetropolisFolder filesep ModelName '_mh*.mat']);
if ~length(mh_files)
error('Estimation::mcmc: I cannot find any MH file to load here!')
end
fidlog = fopen([MetropolisFolder '/metropolis.log'],'a');
fprintf(fidlog,'\n');
fprintf(fidlog,['%% Session ' int2str(length(record.MhDraws(:,1))+1) '.\n']);
fprintf(fidlog,' \n');
fprintf(fidlog,[' Number of blocks...............: ' int2str(nblck) '\n']);
fprintf(fidlog,[' Number of simulations per block: ' int2str(nruns(1)) '\n']);
fprintf(fidlog,' \n');
past_number_of_blocks = record.Nblck;
if past_number_of_blocks ~= nblck
disp('Estimation::mcmc: The specified number of blocks doesn''t match with the previous number of blocks!')
disp(['Estimation::mcmc: You declared ' int2str(nblck) ' blocks, but the previous number of blocks was ' int2str(past_number_of_blocks) '.'])
disp(['Estimation::mcmc: I will run the Metropolis-Hastings with ' int2str(past_number_of_blocks) ' blocks.' ])
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:
LastFileNumber = record.LastFileNumber;
LastLineNumber = record.LastLineNumber;
if LastLineNumber < MAX_nruns
NewFile = ones(nblck,1)*LastFileNumber;
fline = ones(nblck,1)*(LastLineNumber+1);
else
NewFile = ones(nblck,1)*(LastFileNumber+1);
fline = ones(nblck,1);
end
ilogpo2 = record.LastLogPost;
ix2 = record.LastParameters;
fblck = 1;
NumberOfPreviousSimulations = sum(record.MhDraws(:,1),1);
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;
AnticipatedNumberOfFiles = ceil(NumberOfDrawsToBeSaved/MAX_nruns);
AnticipatedNumberOfLinesInTheLastFile = NumberOfDrawsToBeSaved - (AnticipatedNumberOfFiles-1)*MAX_nruns;
record.LastFileNumber = LastFileNumber + AnticipatedNumberOfFiles;
record.LastLineNumber = AnticipatedNumberOfLinesInTheLastFile;
record.MhDraws(end,1) = nruns(1);
record.MhDraws(end,2) = AnticipatedNumberOfFiles;
record.MhDraws(end,3) = AnticipatedNumberOfLinesInTheLastFile;
record.InitialSeeds = record.LastSeeds;
write_mh_history_file(MetropolisFolder, ModelName, record);
fprintf('Done.\n')
disp(['Estimation::mcmc: Ok. I have loaded ' int2str(NumberOfPreviousSimulations) ' simulations.'])
skipline()
fclose(fidlog);
elseif options_.mh_recover
% The previous metropolis-hastings crashed before the end! I try to recover the saved draws...
disp('Estimation::mcmc: Recover mode!')
load_last_mh_history_file(MetropolisFolder, ModelName);
nblck = record.Nblck;% Number of "parallel" mcmc chains.
options_.mh_nblck = nblck;
if size(record.MhDraws,1) == 1
OldMh = 0;% The crashed metropolis was the first session.
else
OldMh = 1;% The crashed metropolis wasn't the first session.
end
% Default initialization:
if OldMh
ilogpo2 = record.LastLogPost;
ix2 = record.LastParameters;
else
ilogpo2 = record.InitialLogPost;
ix2 = record.InitialParameters;
end
% Set NewFile, a nblck*1 vector of integers, and fline (first line), a nblck*1 vector of integers.
if OldMh
LastLineNumberInThePreviousMh = record.MhDraws(end-1,3);% Number of lines in the last mh files of the previous session.
LastFileNumberInThePreviousMh = sum(record.MhDraws(1:end-1,2),1);% Number of mh files in the the previous sessions.
if LastLineNumberInThePreviousMh < MAX_nruns% Test if the last mh files of the previous session are not complete (yes)
NewFile = ones(nblck,1)*LastFileNumberInThePreviousMh;
fline = ones(nblck,1)*(LastLineNumberInThePreviousMh+1);
else% The last mh files of the previous session are complete.
NewFile = ones(nblck,1)*(LastFileNumberInThePreviousMh+1);
fline = ones(nblck,1);
end
else
LastLineNumberInThePreviousMh = 0;
LastFileNumberInThePreviousMh = 0;
NewFile = ones(nblck,1);
fline = ones(nblck,1);
end
% Set fblck (First block), an integer targeting the crashed mcmc chain.
fblck = 1;
% How many mh files should we have ?
ExpectedNumberOfMhFilesPerBlock = sum(record.MhDraws(:,2),1);
ExpectedNumberOfMhFiles = ExpectedNumberOfMhFilesPerBlock*nblck;
% I count the total number of saved mh files...
AllMhFiles = dir([BaseName '_mh*_blck*.mat']);
TotalNumberOfMhFiles = length(AllMhFiles);
% And I quit if I can't find a crashed mcmc chain.
if (TotalNumberOfMhFiles==ExpectedNumberOfMhFiles)
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('Estimation::mcmc: mh_recover option not required!')
end
% I count the number of saved mh files per block.
NumberOfMhFilesPerBlock = zeros(nblck,1);
for b = 1:nblck
BlckMhFiles = dir([BaseName '_mh*_blck' int2str(b) '.mat']);
NumberOfMhFilesPerBlock(b) = length(BlckMhFiles);
end
% Is there a chain with less mh files than expected ?
while fblck <= nblck
if NumberOfMhFilesPerBlock(fblck) < ExpectedNumberOfMhFilesPerBlock
disp(['Estimation::mcmc: Chain ' int2str(fblck) ' is not complete!'])
break
% The mh_recover session will start from chain fblck.
else
disp(['Estimation::mcmc: Chain ' int2str(fblck) ' is complete!'])
end
fblck = fblck+1;
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
% that NumberOfSavedMhFilesInTheCrashedBlck is the number of saved mh files in the crashed chain
% of the current session).
if OldMh
NumberOfSavedMhFilesInTheCrashedBlck = NumberOfSavedMhFilesInTheCrashedBlck - LastFileNumberInThePreviousMh;
end
NumberOfSavedMhFiles = NumberOfSavedMhFilesInTheCrashedBlck+LastFileNumberInThePreviousMh;
% Correct initial conditions.
if NumberOfSavedMhFiles
load([BaseName '_mh' int2str(NumberOfSavedMhFiles) '_blck' int2str(fblck) '.mat']);
ilogpo2(fblck) = logpo2(end);
ix2(fblck,:) = x2(end,:);
end
end

93
matlab/posterior_sampler_iteration.m Normal file
View File

@ -0,0 +1,93 @@
function [par, logpost, accepted, neval] = posterior_sampler_iteration(TargetFun,last_draw, last_posterior, sampler_options,varargin)
% function [par, logpost, accepted, neval] = posterior_sampler_iteration(TargetFun,last_draw, last_posterior, sampler_options,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_)
% posterior samplers
%
% INPUTS
% posterior_sampler_options: posterior sampler options
% options_: structure storing the options
% OUTPUTS
% posterior_sampler_options: checked posterior sampler options
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2015 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 <http://www.gnu.org/licenses/>.
posterior_sampling_method = sampler_options.posterior_sampling_method;
mh_bounds = sampler_options.bounds;
switch posterior_sampling_method
case 'slice'
[par, logpost, neval] = slice_sampler(TargetFun,last_draw, [mh_bounds.lb mh_bounds.ub], sampler_options,varargin{:});
accepted = 1;
case 'random_walk_metropolis_hastings'
neval = 1;
ProposalFun = sampler_options.proposal_distribution;
proposal_covariance_Cholesky_decomposition = sampler_options.proposal_covariance_Cholesky_decomposition;
n = sampler_options.n;
par = feval(ProposalFun, last_draw, proposal_covariance_Cholesky_decomposition, n);
if all( par(:) > mh_bounds.lb ) && all( par(:) < mh_bounds.ub )
try
logpost = - feval(TargetFun, par(:),varargin{:});
catch
logpost = -inf;
end
else
logpost = -inf;
end
r = logpost-last_posterior;
if (logpost > -inf) && (log(rand) < r)
accepted = 1;
else
accepted = 0;
par = last_draw;
logpost = last_posterior;
end
case 'independent_metropolis_hastings'
neval = 1;
ProposalFun = sampler_options.proposal_distribution;
ProposalDensity = sampler_options.ProposalDensity;
proposal_covariance_Cholesky_decomposition = sampler_options.proposal_covariance_Cholesky_decomposition;
n = sampler_options.n;
xparam1 = sampler_options.xparam1;
par = feval(ProposalFun, sampler_options.xparam1, proposal_covariance_Cholesky_decomposition, n);
if all( par(:) > mh_bounds.lb ) && all( par(:) < mh_bounds.ub )
try
logpost = - feval(TargetFun, par(:),varargin{:});
catch
logpost = -inf;
end
else
logpost = -inf;
end
r = logpost - last_posterior + ...
log(feval(ProposalDensity, last_draw, xparam1, proposal_covariance_Cholesky_decomposition, n)) - ...
log(feval(ProposalDensity, par, xparam1, proposal_covariance_Cholesky_decomposition, n));
if (logpost > -inf) && (log(rand) < r)
accepted = 1;
else
accepted = 0;
par = last_draw;
logpost = last_posterior;
end
end