posterior_sampler.m 8.65 KB
Newer Older
1 2
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_)
3 4 5
% Random Walk Metropolis-Hastings algorithm.
%
% INPUTS
6 7 8 9 10 11 12
%   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).
13
%   o mh_bounds  [double]   (p*2) matrix defining lower and upper bounds for the parameters.
14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
%   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
%
% 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)
29
% has been removed from this function and been placed in the posterior_sampler_core.m funtion.
30
%
31 32 33 34 35 36 37 38 39
% 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.

40
% Copyright (C) 2006-2017 Dynare Team
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57
%
% 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/>.

vv = sampler_options.invhess;
58
% Initialization of the sampler
59
[ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d, bayestopt_] = ...
60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109
    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.

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',[]);

110
if strcmp(sampler_options.posterior_sampling_method,'tailored_random_block_metropolis_hastings')
111
    localVars.options_.silent_optimizer=1; %locally set optimizer to silent mode
112 113 114
    if ~isempty(sampler_options.optim_opt)
        localVars.options_.optim_opt=sampler_options.optim_opt; %locally set options for optimizer
    end
115
end
116 117

% User doesn't want to use parallel computing, or wants to compute a
118
% single chain compute sequentially.
119

120
if isnumeric(options_.parallel) || (nblck-fblck)==0
121
    fout = posterior_sampler_core(localVars, fblck, nblck, 0);
122
    record = fout.record;
123 124
    % Parallel in Local or remote machine.
else
125 126 127 128 129
    % 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']};
130 131 132
    if M_.set_auxiliary_variables
        NamFileInput(3,:) = {'',[M_.fname '_set_auxiliary_variables.m']};
    end
133 134
    if options_.steadystate_flag
        if options_.steadystate_flag == 1
135 136 137 138
            NamFileInput(length(NamFileInput)+1,:)={'',[M_.fname '_steadystate.m']};
        else
            NamFileInput(length(NamFileInput)+1,:)={'',[M_.fname '_steadystate2.m']};
        end
139
    end
140
    if (options_.load_mh_file~=0)  && any(fline>1)
141 142 143 144
        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/'],'*.*'};
145
    [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'posterior_sampler_core', localVars, globalVars, options_.parallel_info);
146
    for j=1:totCPU
147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185
        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
186
end