independent_metropolis_hastings_core

 PARALLEL CONTEXT
 The most computationally intensive portion of code in
 independent_metropolis_hastings (the 'for xxx = fblck:nblck' cycle).
 See the comment in random_walk_metropolis_hastings_core.m funtion.

 INPUTS
   See See the comment in random_walk_metropolis_hastings_core.m funtion.

0001 function myoutput = independent_metropolis_hastings_core(myinputs,fblck,nblck,whoiam, ThisMatlab)
0002 % PARALLEL CONTEXT
0003 % The most computationally intensive portion of code in
0004 % independent_metropolis_hastings (the 'for xxx = fblck:nblck' cycle).
0005 % See the comment in random_walk_metropolis_hastings_core.m funtion.
0006 %
0007 % INPUTS
0008 %   See See the comment in random_walk_metropolis_hastings_core.m funtion.
0009 
0010 % OUTPUTS
0011 %   See See the comment in random_walk_metropolis_hastings_core.m funtion.
0012 %
0013 % ALGORITHM
0014 %   Portion of Independing Metropolis-Hastings.
0015 %
0016 % SPECIAL REQUIREMENTS.
0017 %   None.
0018 %
0019 % Copyright (C) 2006-2011 Dynare Team
0020 %
0021 % This file is part of Dynare.
0022 %
0023 % Dynare is free software: you can redistribute it and/or modify
0024 % it under the terms of the GNU General Public License as published by
0025 % the Free Software Foundation, either version 3 of the License, or
0026 % (at your option) any later version.
0027 %
0028 % Dynare is distributed in the hope that it will be useful,
0029 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0030 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0031 % GNU General Public License for more details.
0032 %
0033 % You should have received a copy of the GNU General Public License
0034 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0035 
0036 if nargin<4,
0037     whoiam=0;
0038 end
0039 
0040 
0041 global bayestopt_ estim_params_ options_  M_ oo_
0042 
0043 % Reshape 'myinputs' for local computation.
0044 % In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by:
0045 
0046 TargetFun=myinputs.TargetFun;
0047 ProposalFun=myinputs.ProposalFun;
0048 xparam1=myinputs.xparam1;
0049 vv=myinputs.vv;
0050 mh_bounds=myinputs.mh_bounds;
0051 ix2=myinputs.ix2;
0052 ilogpo2=myinputs.ilogpo2;
0053 ModelName=myinputs.ModelName;
0054 fline=myinputs.fline;
0055 npar=myinputs.npar;
0056 nruns=myinputs.nruns;
0057 NewFile=myinputs.NewFile;
0058 MAX_nruns=myinputs.MAX_nruns;
0059 d=myinputs.d;
0060 InitSizeArray=myinputs.InitSizeArray;
0061 record=myinputs.record;
0062 varargin=myinputs.varargin;
0063 
0064 if whoiam
0065     Parallel=myinputs.Parallel;
0066     % initialize persistent variables in priordens()
0067     priordens(xparam1',bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, ...
0068               bayestopt_.p3,bayestopt_.p4,1);
0069 end
0070 
0071 % (re)Set the penalty.
0072 bayestopt_.penalty = Inf;
0073 
0074 MhDirectoryName = CheckPath('metropolis',M_.dname);
0075 
0076 OpenOldFile = ones(nblck,1);
0077 if strcmpi(ProposalFun,'rand_multivariate_normal')
0078     n = npar;
0079     ProposalDensity = 'multivariate_normal_pdf';
0080 elseif strcmpi(ProposalFun,'rand_multivariate_student')
0081     n = options_.student_degrees_of_freedom;
0082     ProposalDensity = 'multivariate_student_pdf';
0083 end
0084 % load([MhDirectoryName '/' ModelName '_mh_history.mat'],'record');
0085 %%%%
0086 %%%% NOW i run the (nblck-fblck+1) metropolis-hastings chains
0087 %%%%
0088 
0089 if any(isnan(bayestopt_.jscale))
0090     if exist([ModelName '_optimal_mh_scale_parameter.mat'])% This file is created by mode_compute=6.
0091         load([ModelName '_optimal_mh_scale_parameter'])
0092         proposal_covariance = d*Scale;
0093     else
0094         error('mh:: Something is wrong. I can''t figure out the value of the scale parameter.')
0095     end
0096 else
0097     proposal_covariance = d*diag(bayestopt_.jscale);
0098 end
0099 
0100 jloop=0;
0101 
0102 for b = fblck:nblck,
0103     jloop=jloop+1;
0104     randn('state',record.Seeds(b).Normal);
0105     rand('state',record.Seeds(b).Unifor);
0106     if (options_.load_mh_file~=0)  && (fline(b)>1) && OpenOldFile(b)
0107         load(['./' MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) ...
0108               '_blck' int2str(b) '.mat'])
0109         x2 = [x2;zeros(InitSizeArray(b)-fline(b)+1,npar)];
0110         logpo2 = [logpo2;zeros(InitSizeArray(b)-fline(b)+1,1)];
0111         OpenOldFile(b) = 0;
0112     else
0113         x2 = zeros(InitSizeArray(b),npar);
0114         logpo2 = zeros(InitSizeArray(b),1);
0115     end
0116     if exist('OCTAVE_VERSION') || options_.console_mode
0117         diary off
0118         disp(' ')
0119     elseif whoiam
0120         %       keyboard;
0121         waitbarString = ['Please wait... Metropolis-Hastings (' int2str(b) '/' int2str(options_.mh_nblck) ')...'];
0122         %       waitbarTitle=['Metropolis-Hastings ',options_.parallel(ThisMatlab).ComputerName];
0123         if options_.parallel(ThisMatlab).Local,
0124             waitbarTitle=['Local '];
0125         else
0126             waitbarTitle=[options_.parallel(ThisMatlab).ComputerName];
0127         end
0128         fMessageStatus(0,whoiam,waitbarString, waitbarTitle, options_.parallel(ThisMatlab));
0129     else,
0130         hh = waitbar(0,['Please wait... Metropolis-Hastings (' int2str(b) '/' int2str(options_.mh_nblck) ')...']);
0131         set(hh,'Name','Metropolis-Hastings');
0132 
0133     end
0134     isux = 0;
0135     jsux = 0;
0136     irun = fline(b);
0137     j = 1;
0138     while j <= nruns(b)
0139         par = feval(ProposalFun, xparam1, proposal_covariance, n);
0140         if all( par(:) > mh_bounds(:,1) ) && all( par(:) < mh_bounds(:,2) )
0141             try
0142                 logpost = - feval(TargetFun, par(:),varargin{:});
0143             catch,
0144                 logpost = -inf;
0145             end
0146         else
0147             logpost = -inf;
0148         end
0149         r = logpost - ilogpo2(b) + ...
0150             log(feval(ProposalDensity, ix2(b,:), xparam1, proposal_covariance, n)) - ...
0151             log(feval(ProposalDensity, par, xparam1, proposal_covariance, n));
0152         if (logpost > -inf) && (log(rand) < r)
0153             x2(irun,:) = par;
0154             ix2(b,:) = par;
0155             logpo2(irun) = logpost;
0156             ilogpo2(b) = logpost;
0157             isux = isux + 1;
0158             jsux = jsux + 1;
0159         else
0160             x2(irun,:) = ix2(b,:);
0161             logpo2(irun) = ilogpo2(b);
0162         end
0163         prtfrc = j/nruns(b);
0164         if exist('OCTAVE_VERSION') || options_.console_mode
0165             if mod(j, 10) == 0
0166                 if exist('OCTAVE_VERSION')
0167                     if (whoiam==0),
0168                         printf('MH: Computing Metropolis-Hastings (chain %d/%d): %3.f%% done, acception rate: %3.f%%\r', b, nblck, 100 * prtfrc, 100 * isux / j);
0169                     end
0170                 else
0171                     fprintf('   MH: Computing Metropolis-Hastings (chain %d/%d): %3.f \b%% done, acception rate: %3.f \b%%\r', b, nblck, 100 * prtfrc, 100 * isux / j);
0172                 end
0173             end
0174             if mod(j,50)==0 && whoiam,
0175                 %             keyboard;
0176                 waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) '), ' sprintf('accept. %3.f%%%%', 100 * isux/j)];
0177                 fMessageStatus(prtfrc,whoiam,waitbarString, '', options_.parallel(ThisMatlab))
0178             end
0179         else
0180             if mod(j, 3)==0 && ~whoiam
0181                 waitbar(prtfrc,hh,[ '(' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)]);
0182             elseif mod(j,50)==0 && whoiam,
0183                 %             keyboard;
0184                 waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)];
0185                 fMessageStatus(prtfrc,whoiam,waitbarString, waitbarTitle, options_.parallel(ThisMatlab))
0186             end
0187         end
0188 
0189         if (irun == InitSizeArray(b)) || (j == nruns(b)) % Now I save the simulations
0190             save([MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) '_blck' int2str(b) '.mat'],'x2','logpo2');
0191             fidlog = fopen([MhDirectoryName '/metropolis.log'],'a');
0192             fprintf(fidlog,['\n']);
0193             fprintf(fidlog,['%% Mh' int2str(NewFile(b)) 'Blck' int2str(b) ' (' datestr(now,0) ')\n']);
0194             fprintf(fidlog,' \n');
0195             fprintf(fidlog,['  Number of simulations.: ' int2str(length(logpo2)) '\n']);
0196             fprintf(fidlog,['  Acceptation rate......: ' num2str(jsux/length(logpo2)) '\n']);
0197             fprintf(fidlog,['  Posterior mean........:\n']);
0198             for i=1:length(x2(1,:))
0199                 fprintf(fidlog,['    params:' int2str(i) ': ' num2str(mean(x2(:,i))) '\n']);
0200             end
0201             fprintf(fidlog,['    log2po:' num2str(mean(logpo2)) '\n']);
0202             fprintf(fidlog,['  Minimum value.........:\n']);;
0203             for i=1:length(x2(1,:))
0204                 fprintf(fidlog,['    params:' int2str(i) ': ' num2str(min(x2(:,i))) '\n']);
0205             end
0206             fprintf(fidlog,['    log2po:' num2str(min(logpo2)) '\n']);
0207             fprintf(fidlog,['  Maximum value.........:\n']);
0208             for i=1:length(x2(1,:))
0209                 fprintf(fidlog,['    params:' int2str(i) ': ' num2str(max(x2(:,i))) '\n']);
0210             end
0211             fprintf(fidlog,['    log2po:' num2str(max(logpo2)) '\n']);
0212             fprintf(fidlog,' \n');
0213             fclose(fidlog);
0214             jsux = 0;
0215             if j == nruns(b) % I record the last draw...
0216                 record.LastParameters(b,:) = x2(end,:);
0217                 record.LastLogLiK(b) = logpo2(end);
0218             end
0219             % size of next file in chain b
0220             InitSizeArray(b) = min(nruns(b)-j,MAX_nruns);
0221             % initialization of next file if necessary
0222             if InitSizeArray(b)
0223                 x2 = zeros(InitSizeArray(b),npar);
0224                 logpo2 = zeros(InitSizeArray(b),1);
0225                 NewFile(b) = NewFile(b) + 1;
0226                 irun = 0;
0227             end
0228         end
0229         j=j+1;
0230         irun = irun + 1;
0231     end% End of the simulations for one mh-block.
0232     record.AcceptationRates(b) = isux/j;
0233     if exist('OCTAVE_VERSION') || options_.console_mode
0234         printf('\n');
0235         diary on
0236     elseif ~whoiam
0237         close(hh);
0238     end
0239     record.Seeds(b).Normal = randn('state');
0240     record.Seeds(b).Unifor = rand('state');
0241     OutputFileName(jloop,:) = {[MhDirectoryName,filesep], [ModelName '_mh*_blck' int2str(b) '.mat']};
0242 end% End of the loop over the mh-blocks.
0243 
0244 
0245 myoutput.record = record;
0246 myoutput.irun = irun;
0247 myoutput.NewFile = NewFile;
0248 myoutput.OutputFileName = OutputFileName;
0249 
0250