diff --git a/matlab/check_posterior_sampler_options.m b/matlab/check_posterior_sampler_options.m
new file mode 100644
index 0000000000000000000000000000000000000000..5f0cd7b7d39f88527bd81181b17d5c137dc21aac
--- /dev/null
+++ b/matlab/check_posterior_sampler_options.m
@@ -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
+
+
+
diff --git a/matlab/posterior_sampler.m b/matlab/posterior_sampler.m
new file mode 100644
index 0000000000000000000000000000000000000000..4ec8aa14066d769b6dd41f1ccd73813556059482
--- /dev/null
+++ b/matlab/posterior_sampler.m
@@ -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
\ No newline at end of file
diff --git a/matlab/posterior_sampler_core.m b/matlab/posterior_sampler_core.m
new file mode 100644
index 0000000000000000000000000000000000000000..d8a0dec712017b58bd9d6284c8dc8be6a6417ed0
--- /dev/null
+++ b/matlab/posterior_sampler_core.m
@@ -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;
\ No newline at end of file
diff --git a/matlab/posterior_sampler_initialization.m b/matlab/posterior_sampler_initialization.m
new file mode 100644
index 0000000000000000000000000000000000000000..2e3804453468d2c85c1a7e056ad92307ca0323f3
--- /dev/null
+++ b/matlab/posterior_sampler_initialization.m
@@ -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
\ No newline at end of file
diff --git a/matlab/posterior_sampler_iteration.m b/matlab/posterior_sampler_iteration.m
new file mode 100644
index 0000000000000000000000000000000000000000..bbd642d750dfdbb0d87b34cdd5012093654e3b4e
--- /dev/null
+++ b/matlab/posterior_sampler_iteration.m
@@ -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
\ No newline at end of file