independent_metropolis_hastings

 Independent 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 
   None  

 ALGORITHM 
   Metropolis-Hastings.       

 SPECIAL REQUIREMENTS
   None.

 PARALLEL CONTEXT
 See the comment in random_walk_metropolis_hastings.m funtion.

0001 function independent_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin)
0002 
0003 % Independent 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 %   None
0015 %
0016 % ALGORITHM
0017 %   Metropolis-Hastings.
0018 %
0019 % SPECIAL REQUIREMENTS
0020 %   None.
0021 %
0022 % PARALLEL CONTEXT
0023 % See the comment in random_walk_metropolis_hastings.m funtion.
0024 
0025 % Copyright (C) 2006-2011 Dynare Team
0026 %
0027 % This file is part of Dynare.
0028 %
0029 % Dynare is free software: you can redistribute it and/or modify
0030 % it under the terms of the GNU General Public License as published by
0031 % the Free Software Foundation, either version 3 of the License, or
0032 % (at your option) any later version.
0033 %
0034 % Dynare is distributed in the hope that it will be useful,
0035 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0036 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0037 % GNU General Public License for more details.
0038 %
0039 % You should have received a copy of the GNU General Public License
0040 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0041 
0042 global M_ options_ bayestopt_ estim_params_ oo_
0043 %%%%
0044 %%%% Initialization of the independent metropolis-hastings chains.
0045 %%%%
0046 [ ix2, ilogpo2, ModelName, MhDirectoryName, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ...
0047     metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds, varargin{:});
0048 
0049 xparam1 = transpose(xparam1);      
0050 InitSizeArray = min([repmat(MAX_nruns,nblck,1) fline+nruns-1],[],2);
0051 
0052 load([MhDirectoryName '/' ModelName '_mh_history.mat'],'record');
0053 
0054 %The mandatory variables for local/remote parallel computing are stored in localVars struct.
0055 
0056 localVars =   struct('TargetFun', TargetFun, ...
0057                      'ProposalFun', ProposalFun, ...
0058                      'xparam1', xparam1, ...
0059                      'vv', vv, ...
0060                      'mh_bounds', mh_bounds, ...
0061                      'ix2', ix2, ...
0062                      'ilogpo2', ilogpo2, ...
0063                      'ModelName', ModelName, ...
0064                      'fline', fline, ...
0065                      'npar', npar, ...
0066                      'nruns', nruns, ...
0067                      'NewFile', NewFile, ...
0068                      'MAX_nruns', MAX_nruns, ...
0069                      'd', d);
0070 localVars.InitSizeArray=InitSizeArray;
0071 localVars.record=record;
0072 localVars.varargin=varargin;
0073 
0074 % Like a sequential execution!
0075 if isnumeric(options_.parallel),
0076     fout = independent_metropolis_hastings_core(localVars, fblck, nblck, 0);
0077     record = fout.record;
0078     % Parallel execution.
0079 else
0080     % global variables for parallel routines
0081     globalVars = struct('M_',M_, ...
0082                         'options_', options_, ...
0083                         'bayestopt_', bayestopt_, ...
0084                         'estim_params_', estim_params_, ...
0085                         'oo_', oo_);
0086     
0087     % which files have to be copied to run remotely
0088     NamFileInput(1,:) = {'',[ModelName '_static.m']};
0089     NamFileInput(2,:) = {'',[ModelName '_dynamic.m']};
0090     if options_.steadystate_flag,
0091         NamFileInput(length(NamFileInput)+1,:)={'',[ModelName '_steadystate.m']};
0092     end
0093     if (options_.load_mh_file~=0) && any(fline>1) ,
0094         NamFileInput(length(NamFileInput)+1,:)={[M_.dname '/metropolis/'],[ModelName '_mh' int2str(NewFile(1)) '_blck*.mat']};
0095     end
0096     if exist([ModelName '_optimal_mh_scale_parameter.mat'])
0097         NamFileInput(length(NamFileInput)+1,:)={'',[ModelName '_optimal_mh_scale_parameter.mat']};
0098     end
0099     
0100     % from where to get back results
0101     %     NamFileOutput(1,:) = {[M_.dname,'/metropolis/'],'*.*'};
0102     
0103     [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'independent_metropolis_hastings_core', localVars, globalVars, options_.parallel_info);
0104     for j=1:totCPU,
0105         offset = sum(nBlockPerCPU(1:j-1))+fblck-1;
0106         record.LastLogLiK(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastLogLiK(offset+1:sum(nBlockPerCPU(1:j)));
0107         record.LastParameters(offset+1:sum(nBlockPerCPU(1:j)),:)=fout(j).record.LastParameters(offset+1:sum(nBlockPerCPU(1:j)),:);
0108         record.AcceptationRates(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.AcceptationRates(offset+1:sum(nBlockPerCPU(1:j)));
0109         record.Seeds(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.Seeds(offset+1:sum(nBlockPerCPU(1:j)));
0110     end
0111 
0112 end
0113 
0114 irun = fout(1).irun;
0115 NewFile = fout(1).NewFile;
0116 
0117 % record.Seeds.Normal = randn('state');
0118 % record.Seeds.Unifor = rand('state');
0119 save([MhDirectoryName '/' ModelName '_mh_history.mat'],'record');
0120 disp(['MH: Number of mh files                   : ' int2str(NewFile(1)) ' per block.'])
0121 disp(['MH: Total number of generated files      : ' int2str(NewFile(1)*nblck) '.'])
0122 disp(['MH: Total number of iterations           : ' int2str((NewFile(1)-1)*MAX_nruns+irun-1) '.'])
0123 disp('MH: average acceptation rate per chain   : ')
0124 disp(record.AcceptationRates);
0125 disp(' ')