1function posterior_sampler(TargetFun,ProposalFun,xparam1,sampler_options,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_) 2% function posterior_sampler(TargetFun,ProposalFun,xparam1,sampler_options,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_) 3% Random Walk Metropolis-Hastings algorithm. 4% 5% INPUTS 6% o TargetFun [char] string specifying the name of the objective 7% function (posterior kernel). 8% o ProposalFun [char] string specifying the name of the proposal 9% density 10% o xparam1 [double] (p*1) vector of parameters to be estimated (initial values). 11% o sampler_options structure 12% .invhess [double] (p*p) matrix, posterior covariance matrix (at the mode). 13% o mh_bounds [double] (p*2) matrix defining lower and upper bounds for the parameters. 14% o dataset_ data structure 15% o dataset_info dataset info structure 16% o options_ options structure 17% o M_ model structure 18% o estim_params_ estimated parameters structure 19% o bayestopt_ estimation options structure 20% o oo_ outputs structure 21% 22% SPECIAL REQUIREMENTS 23% None. 24% 25% PARALLEL CONTEXT 26% The most computationally intensive part of this function may be executed 27% in parallel. The code suitable to be executed in 28% parallel on multi core or cluster machine (in general a 'for' cycle) 29% has been removed from this function and been placed in the posterior_sampler_core.m funtion. 30% 31% The DYNARE parallel packages comprise a i) set of pairs of Matlab functions that can be executed in 32% parallel and called name_function.m and name_function_core.m and ii) a second set of functions used 33% to manage the parallel computations. 34% 35% This function was the first function to be parallelized. Later, other 36% functions have been parallelized using the same methodology. 37% Then the comments write here can be used for all the other pairs of 38% parallel functions and also for management functions. 39 40% Copyright (C) 2006-2017 Dynare Team 41% 42% This file is part of Dynare. 43% 44% Dynare is free software: you can redistribute it and/or modify 45% it under the terms of the GNU General Public License as published by 46% the Free Software Foundation, either version 3 of the License, or 47% (at your option) any later version. 48% 49% Dynare is distributed in the hope that it will be useful, 50% but WITHOUT ANY WARRANTY; without even the implied warranty of 51% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 52% GNU General Public License for more details. 53% 54% You should have received a copy of the GNU General Public License 55% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 56 57vv = sampler_options.invhess; 58% Initialization of the sampler 59[ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d, bayestopt_] = ... 60 posterior_sampler_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); 61 62InitSizeArray = min([repmat(MAX_nruns,nblck,1) fline+nruns-1],[],2); 63 64% Load last mh history file 65load_last_mh_history_file(MetropolisFolder, ModelName); 66 67% Only for test parallel results!!! 68 69% To check the equivalence between parallel and serial computation! 70% First run in serial mode, and then comment the follow line. 71% save('recordSerial.mat','-struct', 'record'); 72 73% For parallel runs after serial runs with the abobe line active. 74% TempRecord=load('recordSerial.mat'); 75% record.Seeds=TempRecord.Seeds; 76 77 78 79% Snapshot of the current state of computing. It necessary for the parallel 80% execution (i.e. to execute in a corretct way a portion of code remotely or 81% on many cores). The mandatory variables for local/remote parallel 82% computing are stored in the localVars struct. 83 84localVars = struct('TargetFun', TargetFun, ... 85 'ProposalFun', ProposalFun, ... 86 'xparam1', xparam1, ... 87 'vv', vv, ... 88 'sampler_options', sampler_options, ... 89 'mh_bounds', mh_bounds, ... 90 'ix2', ix2, ... 91 'ilogpo2', ilogpo2, ... 92 'ModelName', ModelName, ... 93 'fline', fline, ... 94 'npar', npar, ... 95 'nruns', nruns, ... 96 'NewFile', NewFile, ... 97 'MAX_nruns', MAX_nruns, ... 98 'd', d, ... 99 'InitSizeArray',InitSizeArray, ... 100 'record', record, ... 101 'dataset_', dataset_, ... 102 'dataset_info', dataset_info, ... 103 'options_', options_, ... 104 'M_',M_, ... 105 'bayestopt_', bayestopt_, ... 106 'estim_params_', estim_params_, ... 107 'oo_', oo_,... 108 'varargin',[]); 109 110if strcmp(sampler_options.posterior_sampling_method,'tailored_random_block_metropolis_hastings') 111 localVars.options_.silent_optimizer=1; %locally set optimizer to silent mode 112 if ~isempty(sampler_options.optim_opt) 113 localVars.options_.optim_opt=sampler_options.optim_opt; %locally set options for optimizer 114 end 115end 116 117% User doesn't want to use parallel computing, or wants to compute a 118% single chain compute sequentially. 119 120if isnumeric(options_.parallel) || (~isempty(fblck) && (nblck-fblck)==0) 121 fout = posterior_sampler_core(localVars, fblck, nblck, 0); 122 record = fout.record; 123 % Parallel in Local or remote machine. 124else 125 % Global variables for parallel routines. 126 globalVars = struct(); 127 % which files have to be copied to run remotely 128 NamFileInput(1,:) = {'',[ModelName '.static.m']}; 129 NamFileInput(2,:) = {'',[ModelName '.dynamic.m']}; 130 if M_.set_auxiliary_variables 131 NamFileInput(3,:) = {'',[M_.fname '.set_auxiliary_variables.m']}; 132 end 133 if options_.steadystate_flag 134 if options_.steadystate_flag == 1 135 NamFileInput(length(NamFileInput)+1,:)={'',[M_.fname '_steadystate.m']}; 136 else 137 NamFileInput(length(NamFileInput)+1,:)={'',[M_.fname '.steadystate.m']}; 138 end 139 end 140 if (options_.load_mh_file~=0) && any(fline>1) 141 NamFileInput(length(NamFileInput)+1,:)={[M_.dname '/metropolis/'],[ModelName '_mh' int2str(NewFile(1)) '_blck*.mat']}; 142 end 143 % from where to get back results 144 % NamFileOutput(1,:) = {[M_.dname,'/metropolis/'],'*.*'}; 145 if options_.mh_recover && isempty(fblck) 146 % here we just need to retrieve the output of the completed remote jobs 147 fblck=1; 148 options_.parallel_info.parallel_recover = 1; 149 end 150 [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'posterior_sampler_core', localVars, globalVars, options_.parallel_info); 151 for j=1:totCPU 152 offset = sum(nBlockPerCPU(1:j-1))+fblck-1; 153 record.LastLogPost(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastLogPost(offset+1:sum(nBlockPerCPU(1:j))); 154 record.LastParameters(offset+1:sum(nBlockPerCPU(1:j)),:)=fout(j).record.LastParameters(offset+1:sum(nBlockPerCPU(1:j)),:); 155 record.AcceptanceRatio(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.AcceptanceRatio(offset+1:sum(nBlockPerCPU(1:j))); 156 record.FunctionEvalPerIteration(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.FunctionEvalPerIteration(offset+1:sum(nBlockPerCPU(1:j))); 157 record.LastSeeds(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastSeeds(offset+1:sum(nBlockPerCPU(1:j))); 158 end 159 options_.parallel_info.parallel_recover = 0; 160end 161 162irun = fout(1).irun; 163NewFile = fout(1).NewFile; 164 165record.MCMCConcludedSuccessfully = 1; %set indicator for successful run 166 167update_last_mh_history_file(MetropolisFolder, ModelName, record); 168 169% Provide diagnostic output 170skipline() 171disp(['Estimation::mcmc: Number of mh files: ' int2str(NewFile(1)) ' per block.']) 172disp(['Estimation::mcmc: Total number of generated files: ' int2str(NewFile(1)*nblck) '.']) 173disp(['Estimation::mcmc: Total number of iterations: ' int2str((NewFile(1)-1)*MAX_nruns+irun-1) '.']) 174disp(['Estimation::mcmc: Current acceptance ratio per chain: ']) 175for i=1:nblck 176 if i<10 177 disp([' Chain ' num2str(i) ': ' num2str(100*record.AcceptanceRatio(i)) '%']) 178 else 179 disp([' Chain ' num2str(i) ': ' num2str(100*record.AcceptanceRatio(i)) '%']) 180 end 181end 182if max(record.FunctionEvalPerIteration)>1 183 disp(['Estimation::mcmc: Current function evaluations per iteration: ']) 184 for i=1:nblck 185 if i<10 186 disp([' Chain ' num2str(i) ': ' num2str(record.FunctionEvalPerIteration(i))]) 187 else 188 disp([' Chain ' num2str(i) ': ' num2str(record.FunctionEvalPerIteration(i))]) 189 end 190 end 191end 192