2023-09-13 18:09:38 +02:00
function posterior_sampler ( TargetFun,ProposalFun,xparam1,sampler_options,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_,dispString)
% function posterior_sampler(TargetFun,ProposalFun,xparam1,sampler_options,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_,dispString)
2017-05-16 15:10:20 +02:00
% Random Walk Metropolis-Hastings algorithm.
%
% INPUTS
2023-07-13 20:05:52 +02:00
% 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
% o mh_bounds [double] (p*2) matrix defining lower and upper bounds for the parameters.
% o dataset_ [structure] data structure
% o dataset_info [structure] dataset info structure
% o options_ [structure] options structure
% o M_ [structure] model structure
% o estim_params_ [structure] estimated parameters structure
% o bayestopt_ [structure] prior specification structure
% o oo_ [structure] output structure
2023-09-13 18:09:38 +02:00
% o dispString [string] string prependening the messages printed to the command window
2015-09-25 12:54:52 +02:00
%
% 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)
2016-05-19 17:13:39 +02:00
% has been removed from this function and been placed in the posterior_sampler_core.m funtion.
2017-05-16 15:10:20 +02:00
%
2015-09-25 12:54:52 +02:00
% 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.
2023-07-13 18:01:02 +02:00
% Copyright © 2006-2023 Dynare Team
2015-09-25 12:54:52 +02:00
%
% 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
2021-06-09 17:33:48 +02:00
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
2015-09-25 12:54:52 +02:00
vv = sampler_options . invhess ;
2016-05-19 17:13:39 +02:00
% Initialization of the sampler
2016-06-14 15:09:05 +02:00
[ ix2 , ilogpo2 , ModelName , MetropolisFolder , fblck , fline , npar , nblck , nruns , NewFile , MAX_nruns , d , bayestopt_ ] = ...
2023-09-13 18:09:38 +02:00
posterior_sampler_initialization ( TargetFun , xparam1 , vv , mh_bounds , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ , dispString ) ;
2015-09-25 12:54:52 +02:00
InitSizeArray = min ( [ repmat ( MAX_nruns , nblck , 1 ) fline + nruns - 1 ] , [ ] , 2 ) ;
% Load last mh history file
2023-07-13 18:01:02 +02:00
record = load_last_mh_history_file ( MetropolisFolder , ModelName ) ;
2015-09-25 12:54:52 +02:00
% 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.
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' , [ ] ) ;
2017-05-16 12:42:01 +02:00
if strcmp ( sampler_options . posterior_sampling_method , ' tailored_random_block_metropolis_hastings' )
2016-05-13 21:35:59 +02:00
localVars . options_ . silent_optimizer = 1 ; %locally set optimizer to silent mode
2016-05-15 14:18:28 +02:00
if ~ isempty ( sampler_options . optim_opt )
localVars . options_ . optim_opt = sampler_options . optim_opt ; %locally set options for optimizer
end
2016-05-13 21:35:59 +02:00
end
2015-09-25 12:54:52 +02:00
% User doesn't want to use parallel computing, or wants to compute a
2016-05-19 17:13:39 +02:00
% single chain compute sequentially.
2015-09-25 12:54:52 +02:00
2019-07-24 16:16:03 +02:00
if isnumeric ( options_ . parallel ) || ( ~ isempty ( fblck ) && ( nblck - fblck ) == 0 )
2016-04-27 12:14:40 +02:00
fout = posterior_sampler_core ( localVars , fblck , nblck , 0 ) ;
2015-09-25 12:54:52 +02:00
record = fout . record ;
2017-05-16 15:10:20 +02:00
% Parallel in Local or remote machine.
else
2015-09-25 12:54:52 +02:00
% Global variables for parallel routines.
globalVars = struct ( ) ;
% which files have to be copied to run remotely
2018-06-27 17:02:13 +02:00
NamFileInput ( 1 , : ) = { ' ' , [ ModelName ' .static.m' ] } ;
NamFileInput ( 2 , : ) = { ' ' , [ ModelName ' .dynamic.m' ] } ;
2017-12-01 08:48:20 +01:00
if M_ . set_auxiliary_variables
2018-06-27 17:02:13 +02:00
NamFileInput ( 3 , : ) = { ' ' , [ M_ . fname ' .set_auxiliary_variables.m' ] } ;
2017-08-29 10:58:39 +02:00
end
2017-05-16 12:42:01 +02:00
if options_ . steadystate_flag
if options_ . steadystate_flag == 1
2017-02-09 09:44:28 +01:00
NamFileInput ( length ( NamFileInput ) + 1 , : ) = { ' ' , [ M_ . fname ' _steadystate.m' ] } ;
else
2018-06-27 17:02:13 +02:00
NamFileInput ( length ( NamFileInput ) + 1 , : ) = { ' ' , [ M_ . fname ' .steadystate.m' ] } ;
2017-02-09 09:44:28 +01:00
end
2015-09-25 12:54:52 +02:00
end
2017-05-16 12:42:01 +02:00
if ( options_ . load_mh_file ~= 0 ) && any ( fline > 1 )
2015-09-25 12:54:52 +02:00
NamFileInput ( length ( NamFileInput ) + 1 , : ) = { [ M_ . dname ' /metropolis/' ] , [ ModelName ' _mh' int2str ( NewFile ( 1 ) ) ' _blck*.mat' ] } ;
end
% from where to get back results
% NamFileOutput(1,:) = {[M_.dname,'/metropolis/'],'*.*'};
2019-07-24 16:16:03 +02:00
if options_ . mh_recover && isempty ( fblck )
% here we just need to retrieve the output of the completed remote jobs
fblck = 1 ;
options_ . parallel_info . parallel_recover = 1 ;
end
2016-05-13 21:35:59 +02:00
[ fout , nBlockPerCPU , totCPU ] = masterParallel ( options_ . parallel , fblck , nblck , NamFileInput , ' posterior_sampler_core' , localVars , globalVars , options_ . parallel_info ) ;
2017-05-16 12:42:01 +02:00
for j = 1 : totCPU
2015-09-25 12:54:52 +02:00
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 ) ) ) ;
2021-03-25 17:08:18 +01:00
if j == 1
2021-06-08 21:03:08 +02:00
if isfield ( fout ( j ) . record , ' ProposalCovariance' ) && isfield ( fout ( j ) . record , ' ProposalScaleVec' )
record . ProposalCovariance = fout ( j ) . record . ProposalCovariance ;
record . ProposalScaleVec = fout ( j ) . record . ProposalScaleVec ;
end
2021-03-25 17:08:18 +01:00
end
2015-09-25 12:54:52 +02:00
end
2019-07-24 16:16:03 +02:00
options_ . parallel_info . parallel_recover = 0 ;
2015-09-25 12:54:52 +02:00
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 ( )
2023-09-01 21:42:51 +02:00
fprintf ( ' %s: Number of mh files: %u per block.\n' , dispString , NewFile ( 1 ) ) ;
fprintf ( ' %s: Total number of generated files: %u.\n' , dispString , NewFile ( 1 ) * nblck ) ;
fprintf ( ' %s: Total number of iterations: %u.\n' , dispString , ( NewFile ( 1 ) - 1 ) * MAX_nruns + irun - 1 ) ;
fprintf ( ' %s: Current acceptance ratio per chain:\n' , dispString ) ;
2015-09-25 12:54:52 +02:00
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
2023-09-01 21:42:51 +02:00
fprintf ( ' %s: Current function evaluations per iteration:\n' , dispString ) ;
2015-09-25 12:54:52 +02:00
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
2016-05-13 21:35:59 +02:00
end