function record=random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin)
Random walk Metropolis-Hastings algorithm.
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 varargin list of argument following mh_bounds
OUTPUTS
o record [struct] structure describing the iterations
ALGORITHM
Metropolis-Hastings.
SPECIAL REQUIREMENTS
None.
PARALLEL CONTEXT
The most computationally intensive part of this function may be executed
in parallel. The code sutable to be executed in
parallel on multi core or cluster machine (in general a 'for' cycle),
is removed from this function and placed in random_walk_metropolis_hastings_core.m funtion.
Then the DYNARE parallel package contain a set of pairs matlab functions that can be executed in
parallel and called name_function.m and name_function_core.m.
In addition in parallel package we have second set of functions used
to manage the parallel computation.
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 funtions.0001 function record=random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin) 0002 %function record=random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin) 0003 % Random walk Metropolis-Hastings algorithm. 0004 % 0005 % INPUTS 0006 % o TargetFun [char] string specifying the name of the objective 0007 % function (posterior kernel). 0008 % o xparam1 [double] (p*1) vector of parameters to be estimated (initial values). 0009 % o vv [double] (p*p) matrix, posterior covariance matrix (at the mode). 0010 % o mh_bounds [double] (p*2) matrix defining lower and upper bounds for the parameters. 0011 % o varargin list of argument following mh_bounds 0012 % 0013 % OUTPUTS 0014 % o record [struct] structure describing the iterations 0015 % 0016 % ALGORITHM 0017 % Metropolis-Hastings. 0018 % 0019 % SPECIAL REQUIREMENTS 0020 % None. 0021 % 0022 % PARALLEL CONTEXT 0023 % The most computationally intensive part of this function may be executed 0024 % in parallel. The code sutable to be executed in 0025 % parallel on multi core or cluster machine (in general a 'for' cycle), 0026 % is removed from this function and placed in random_walk_metropolis_hastings_core.m funtion. 0027 % Then the DYNARE parallel package contain a set of pairs matlab functions that can be executed in 0028 % parallel and called name_function.m and name_function_core.m. 0029 % In addition in parallel package we have second set of functions used 0030 % to manage the parallel computation. 0031 % 0032 % This function was the first function to be parallelized, later other 0033 % functions have been parallelized using the same methodology. 0034 % Then the comments write here can be used for all the other pairs of 0035 % parallel functions and also for management funtions. 0036 0037 % Copyright (C) 2006-2011 Dynare Team 0038 % 0039 % This file is part of Dynare. 0040 % 0041 % Dynare is free software: you can redistribute it and/or modify 0042 % it under the terms of the GNU General Public License as published by 0043 % the Free Software Foundation, either version 3 of the License, or 0044 % (at your option) any later version. 0045 % 0046 % Dynare is distributed in the hope that it will be useful, 0047 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0048 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0049 % GNU General Public License for more details. 0050 % 0051 % You should have received a copy of the GNU General Public License 0052 % along with Dynare. If not, see <http://www.gnu.org/licenses/>. 0053 0054 global M_ options_ bayestopt_ estim_params_ oo_ 0055 %%%% 0056 %%%% Initialization of the random walk metropolis-hastings chains. 0057 %%%% 0058 [ ix2, ilogpo2, ModelName, MhDirectoryName, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ... 0059 metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds, varargin{:}); 0060 0061 InitSizeArray = min([repmat(MAX_nruns,nblck,1) fline+nruns-1],[],2); 0062 0063 load([MhDirectoryName '/' ModelName '_mh_history.mat'],'record'); 0064 0065 0066 % Only for test parallel results!!! 0067 0068 % To check the equivalence between parallel and serial computation! 0069 % First run in serial mode, and then comment the follow line. 0070 % save('recordSerial.mat','-struct', 'record'); 0071 0072 % For parllel runs after serial runs with the abobe line active. 0073 % TempRecord=load('recordSerial.mat'); 0074 % record.Seeds=TempRecord.Seeds; 0075 0076 0077 0078 % Snapshot of the current state of computing. It necessary for the parallel 0079 % execution (i.e. to execute in a corretct way portion of code remotely or 0080 % on many core). The mandatory variables for local/remote parallel 0081 % computing are stored in localVars struct. 0082 0083 localVars = struct('TargetFun', TargetFun, ... 0084 'ProposalFun', ProposalFun, ... 0085 'xparam1', xparam1, ... 0086 'vv', vv, ... 0087 'mh_bounds', mh_bounds, ... 0088 'ix2', ix2, ... 0089 'ilogpo2', ilogpo2, ... 0090 'ModelName', ModelName, ... 0091 'fline', fline, ... 0092 'npar', npar, ... 0093 'nruns', nruns, ... 0094 'NewFile', NewFile, ... 0095 'MAX_nruns', MAX_nruns, ... 0096 'd', d); 0097 localVars.InitSizeArray=InitSizeArray; 0098 localVars.record=record; 0099 localVars.varargin=varargin; 0100 0101 0102 % The user don't want to use parallel computing, or want to compute a 0103 % single chain. In this cases Random walk Metropolis-Hastings algorithm is 0104 % computed sequentially. 0105 0106 if isnumeric(options_.parallel) || (nblck-fblck)==0, 0107 fout = random_walk_metropolis_hastings_core(localVars, fblck, nblck, 0); 0108 record = fout.record; 0109 0110 % Parallel in Local or remote machine. 0111 else 0112 % Global variables for parallel routines. 0113 globalVars = struct('M_',M_, ... 0114 'options_', options_, ... 0115 'bayestopt_', bayestopt_, ... 0116 'estim_params_', estim_params_, ... 0117 'oo_', oo_); 0118 0119 % which files have to be copied to run remotely 0120 NamFileInput(1,:) = {'',[ModelName '_static.m']}; 0121 NamFileInput(2,:) = {'',[ModelName '_dynamic.m']}; 0122 if options_.steadystate_flag, 0123 NamFileInput(length(NamFileInput)+1,:)={'',[ModelName '_steadystate.m']}; 0124 end 0125 if (options_.load_mh_file~=0) && any(fline>1) , 0126 NamFileInput(length(NamFileInput)+1,:)={[M_.dname '/metropolis/'],[ModelName '_mh' int2str(NewFile(1)) '_blck*.mat']}; 0127 end 0128 if exist([ModelName '_optimal_mh_scale_parameter.mat']) 0129 NamFileInput(length(NamFileInput)+1,:)={'',[ModelName '_optimal_mh_scale_parameter.mat']}; 0130 end 0131 0132 % from where to get back results 0133 % NamFileOutput(1,:) = {[M_.dname,'/metropolis/'],'*.*'}; 0134 0135 [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'random_walk_metropolis_hastings_core', localVars, globalVars, options_.parallel_info); 0136 for j=1:totCPU, 0137 offset = sum(nBlockPerCPU(1:j-1))+fblck-1; 0138 record.LastLogLiK(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastLogLiK(offset+1:sum(nBlockPerCPU(1:j))); 0139 record.LastParameters(offset+1:sum(nBlockPerCPU(1:j)),:)=fout(j).record.LastParameters(offset+1:sum(nBlockPerCPU(1:j)),:); 0140 record.AcceptationRates(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.AcceptationRates(offset+1:sum(nBlockPerCPU(1:j))); 0141 record.Seeds(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.Seeds(offset+1:sum(nBlockPerCPU(1:j))); 0142 end 0143 0144 end 0145 0146 irun = fout(1).irun; 0147 NewFile = fout(1).NewFile; 0148 0149 0150 % record.Seeds.Normal = randn('state'); 0151 % record.Seeds.Unifor = rand('state'); 0152 save([MhDirectoryName '/' ModelName '_mh_history.mat'],'record'); 0153 disp(['MH: Number of mh files : ' int2str(NewFile(1)) ' per block.']) 0154 disp(['MH: Total number of generated files : ' int2str(NewFile(1)*nblck) '.']) 0155 disp(['MH: Total number of iterations : ' int2str((NewFile(1)-1)*MAX_nruns+irun-1) '.']) 0156 disp('MH: average acceptation rate per chain : ') 0157 disp(record.AcceptationRates); 0158 disp(' ')