diff --git a/matlab/independent_metropolis_hastings.m b/matlab/independent_metropolis_hastings.m
index a39b2dc5e..7153c923c 100644
--- a/matlab/independent_metropolis_hastings.m
+++ b/matlab/independent_metropolis_hastings.m
@@ -1,5 +1,6 @@
function independent_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,varargin)
-% Independent Metropolis-Hastings algorithm.
+
+% Independent Metropolis-Hastings algorithm.
%
% INPUTS
% o TargetFun [char] string specifying the name of the objective
@@ -10,10 +11,10 @@ function independent_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bou
% o varargin list of argument following mh_bounds
%
% OUTPUTS
-% None
+% None
%
-% ALGORITHM
-% Metropolis-Hastings.
+% ALGORITHM
+% Metropolis-Hastings.
%
% SPECIAL REQUIREMENTS
% None.
@@ -35,119 +36,88 @@ function independent_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bou
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see .
-global M_ options_ bayestopt_
+global M_ options_ bayestopt_ estim_params_ oo_
%%%%
%%%% Initialization of the independent metropolis-hastings chains.
%%%%
[ ix2, ilogpo2, ModelName, MhDirectoryName, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ...
- metropolis_hastings_initialization(TargetFun,xparam1,vv,mh_bounds,varargin{:});
-xparam1 = transpose(xparam1);
-OpenOldFile = ones(nblck,1);
-if strcmpi(ProposalFun,'rand_multivariate_normal')
- n = npar;
- ProposalDensity = 'multivariate_normal_pdf';
-elseif strcmpi(ProposalFun,'rand_multivariate_student')
- n = options_.student_degrees_of_freedom;
- ProposalDensity = 'multivariate_student_pdf';
-end
-load([MhDirectoryName '/' ModelName '_mh_history'],'record');
-%%%%
-%%%% NOW i run the (nblck-fblck+1) metropolis-hastings chains
-%%%%
-InitSizeArray = min([MAX_nruns*ones(nblck) nruns],[],2);
-jscale = diag(bayestopt_.jscale);
-for b = fblck:nblck
- randn('state',record.Seeds(b).Normal);
- rand('state',record.Seeds(b).Unifor);
- if (options_.load_mh_file~=0) & (fline(b)>1) & OpenOldFile(b)
- load(['./' MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) ...
- '_blck' int2str(b) '.mat'])
- x2 = [x2;zeros(InitSizeArray(b)-fline(b)+1,npar)];
- logpo2 = [logpo2;zeros(InitSizeArray(b)-fline(b)+1,1)];
- OpenOldFile(b) = 0;
- else
- x2 = zeros(InitSizeArray(b),npar);
- logpo2 = zeros(InitSizeArray(b),1);
+ metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds, varargin{:});
+
+xparam1 = transpose(xparam1);
+InitSizeArray = min([repmat(MAX_nruns,nblck,1) fline+nruns-1],[],2);
+
+load([MhDirectoryName '/' ModelName '_mh_history.mat'],'record');
+
+
+localVars = struct('TargetFun', TargetFun, ...
+ 'ProposalFun', ProposalFun, ...
+ 'xparam1', xparam1, ...
+ 'vv', vv, ...
+ 'mh_bounds', mh_bounds, ...
+ 'ix2', ix2, ...
+ 'ilogpo2', ilogpo2, ...
+ 'ModelName', ModelName, ...
+ 'fline', fline, ...
+ 'npar', npar, ...
+ 'nruns', nruns, ...
+ 'NewFile', NewFile, ...
+ 'MAX_nruns', MAX_nruns, ...
+ 'd', d);
+localVars.InitSizeArray=InitSizeArray;
+localVars.record=record;
+localVars.varargin=varargin;
+
+% tic,
+
+
+if isnumeric(options_.parallel),% | isunix, % for the moment exclude unix platform from parallel implementation
+ fout = independent_metropolis_hastings_core(localVars, fblck, nblck, 0);
+ record = fout.record;
+
+else
+ % global variables for parallel routines
+ globalVars = struct('M_',M_, ...
+ 'options_', options_, ...
+ 'bayestopt_', bayestopt_, ...
+ 'estim_params_', estim_params_, ...
+ 'oo_', oo_);
+
+ % which files have to be copied to run remotely
+ NamFileInput(1,:) = {'',[ModelName '_static.m']};
+ NamFileInput(2,:) = {'',[ModelName '_dynamic.m']};
+ if options_.steadystate_flag,
+ NamFileInput(length(NamFileInput)+1,:)={'',[ModelName '_steadystate.m']};
end
- hh = waitbar(0,['Please wait... Metropolis-Hastings (' int2str(b) '/' int2str(nblck) ')...']);
- set(hh,'Name','Metropolis-Hastings');
- isux = 0;
- jsux = 0;
- irun = fline(b);
- j = 1;
- while j <= nruns(b)
- par = feval(ProposalFun, xparam1, d * jscale, n);
- if all(par(:)>mh_bounds(:,1)) && all(par(:) -inf) && (log(rand) < r)
- x2(irun,:) = par;
- ix2(b,:) = par;
- logpo2(irun) = logpost;
- ilogpo2(b) = logpost;
- isux = isux + 1;
- jsux = jsux + 1;
- else
- x2(irun,:) = ix2(b,:);
- logpo2(irun) = ilogpo2(b);
- end
- prtfrc = j/nruns(b);
- waitbar(prtfrc,hh,[ '(' int2str(b) '/' int2str(nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)]);
- if (irun == InitSizeArray(b)) | (j == nruns(b)) % Now I save the simulations
- save([MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) '_blck' int2str(b)],'x2','logpo2');
- InitSizeArray(b) = min(nruns(b)-j,MAX_nruns);
- fidlog = fopen([MhDirectoryName '/metropolis.log'],'a');
- fprintf(fidlog,['\n']);
- fprintf(fidlog,['%% Mh' int2str(NewFile(b)) 'Blck' int2str(b) ' (' datestr(now,0) ')\n']);
- fprintf(fidlog,' \n');
- fprintf(fidlog,[' Number of simulations.: ' int2str(length(logpo2)) '\n']);
- fprintf(fidlog,[' Acceptation rate......: ' num2str(jsux/length(logpo2)) '\n']);
- fprintf(fidlog,[' Posterior mean........:\n']);
- for i=1:length(x2(1,:))
- fprintf(fidlog,[' params:' int2str(i) ': ' num2str(mean(x2(:,i))) '\n']);
- end
- fprintf(fidlog,[' log2po:' num2str(mean(logpo2)) '\n']);
- fprintf(fidlog,[' Minimum value.........:\n']);;
- for i=1:length(x2(1,:))
- fprintf(fidlog,[' params:' int2str(i) ': ' num2str(min(x2(:,i))) '\n']);
- end
- fprintf(fidlog,[' log2po:' num2str(min(logpo2)) '\n']);
- fprintf(fidlog,[' Maximum value.........:\n']);
- for i=1:length(x2(1,:))
- fprintf(fidlog,[' params:' int2str(i) ': ' num2str(max(x2(:,i))) '\n']);
- end
- fprintf(fidlog,[' log2po:' num2str(max(logpo2)) '\n']);
- fprintf(fidlog,' \n');
- fclose(fidlog);
- jsux = 0;
- if j == nruns(b) % I record the last draw...
- record.LastParameters(b,:) = x2(end,:);
- record.LastLogLiK(b) = logpo2(end);
- end
- if InitSizeArray(b)
- x2 = zeros(InitSizeArray(b),npar);
- logpo2 = zeros(InitSizeArray(b),1);
- NewFile(b) = NewFile(b) + 1;
- irun = 0;
- else% InitSizeArray is equal to zero because we are at the end of an mc chain.
- InitSizeArray(b) = min(nruns(b),MAX_nruns);
- end
- end
- j=j+1;
- irun = irun + 1;
- end% End of the simulations for one mh-block.
- record.AcceptationRates(b) = isux/j;
- close(hh);
- record.Seeds(b).Normal = randn('state');
- record.Seeds(b).Unifor = rand('state');
-end% End of the loop over the mh-blocks.
-save([MhDirectoryName '/' ModelName '_mh_history'],'record');
+ if (options_.load_mh_file~=0) & any(fline>1) ,
+ 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/'],'*.*'};
+
+ [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'independent_metropolis_hastings_core', localVars, globalVars);
+ for j=1:totCPU,
+ offset = sum(nBlockPerCPU(1:j-1))+fblck-1;
+ record.LastLogLiK(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastLogLiK(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.AcceptationRates(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.AcceptationRates(offset+1:sum(nBlockPerCPU(1:j)));
+ record.Seeds(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.Seeds(offset+1:sum(nBlockPerCPU(1:j)));
+ end
+
+end
+
+irun = fout(1).irun;
+NewFile = fout(1).NewFile;
+
+
+% ComptationalTime=toc,
+
+% record.Seeds.Normal = randn('state');
+% record.Seeds.Unifor = rand('state');
+save([MhDirectoryName '/' ModelName '_mh_history.mat'],'record');
disp(['MH: Number of mh files : ' int2str(NewFile(1)) ' per block.'])
disp(['MH: Total number of generated files : ' int2str(NewFile(1)*nblck) '.'])
disp(['MH: Total number of iterations : ' int2str((NewFile(1)-1)*MAX_nruns+irun-1) '.'])
+disp('MH: average acceptation rate per chain : ')
+disp(record.AcceptationRates);
disp(' ')
\ No newline at end of file
diff --git a/matlab/independent_metropolis_hastings_core.m b/matlab/independent_metropolis_hastings_core.m
new file mode 100644
index 000000000..a1ec1c05d
--- /dev/null
+++ b/matlab/independent_metropolis_hastings_core.m
@@ -0,0 +1,185 @@
+function myoutput = independent_metropolis_hastings_core(myinputs,fblck,nblck,whoiam, ThisMatlab)
+
+% Copyright (C) 2006-2008 Dynare Team
+%
+% 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 .
+
+if nargin<4,
+ whoiam=0;
+end
+
+
+global bayestopt_ estim_params_ options_ M_ oo_
+
+struct2local(myinputs);
+
+
+MhDirectoryName = CheckPath('metropolis');
+
+OpenOldFile = ones(nblck,1);
+if strcmpi(ProposalFun,'rand_multivariate_normal')
+ n = npar;
+ ProposalDensity = 'multivariate_normal_pdf';
+elseif strcmpi(ProposalFun,'rand_multivariate_student')
+ n = options_.student_degrees_of_freedom;
+ ProposalDensity = 'multivariate_student_pdf';
+end
+% load([MhDirectoryName '/' ModelName '_mh_history.mat'],'record');
+%%%%
+%%%% NOW i run the (nblck-fblck+1) metropolis-hastings chains
+%%%%
+jscale = diag(bayestopt_.jscale);
+
+jloop=0;
+
+for b = fblck:nblck,
+ jloop=jloop+1;
+ randn('state',record.Seeds(b).Normal);
+ rand('state',record.Seeds(b).Unifor);
+ if (options_.load_mh_file~=0) & (fline(b)>1) & OpenOldFile(b)
+ load(['./' MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) ...
+ '_blck' int2str(b) '.mat'])
+ x2 = [x2;zeros(InitSizeArray(b)-fline(b)+1,npar)];
+ logpo2 = [logpo2;zeros(InitSizeArray(b)-fline(b)+1,1)];
+ OpenOldFile(b) = 0;
+ else
+ x2 = zeros(InitSizeArray(b),npar);
+ logpo2 = zeros(InitSizeArray(b),1);
+ end
+ if exist('OCTAVE_VERSION')
+ diary off;
+ elseif whoiam
+% keyboard;
+ waitbarString = ['Please wait... Metropolis-Hastings (' int2str(b) '/' int2str(options_.mh_nblck) ')...'];
+% waitbarTitle=['Metropolis-Hastings ',options_.parallel(ThisMatlab).PcName];
+ if options_.parallel(ThisMatlab).Local,
+ waitbarTitle=['Local '];
+ else
+ waitbarTitle=[options_.parallel(ThisMatlab).PcName];
+ end
+ fMessageStatus(0,whoiam,waitbarString, waitbarTitle, options_.parallel(ThisMatlab), MasterName, DyMo);
+ else,
+ hh = waitbar(0,['Please wait... Metropolis-Hastings (' int2str(b) '/' int2str(options_.mh_nblck) ')...']);
+ set(hh,'Name','Metropolis-Hastings');
+
+ end
+ isux = 0;
+ jsux = 0;
+ irun = fline(b);
+ j = 1;
+ while j <= nruns(b)
+ par = feval(ProposalFun, xparam1, d * jscale, n);
+ if all( par(:) > mh_bounds(:,1) ) & all( par(:) < mh_bounds(:,2) )
+ logpost = - feval(TargetFun, par(:),varargin{:});
+ else
+ logpost = -inf;
+ end
+ r = logpost - ilogpo2(b) + ...
+ log(feval(ProposalDensity, ix2(b,:), xparam1, d, n)) - ...
+ log(feval(ProposalDensity, par, xparam1, d, n));
+ if (logpost > -inf) && (log(rand) < r)
+ x2(irun,:) = par;
+ ix2(b,:) = par;
+ logpo2(irun) = logpost;
+ ilogpo2(b) = logpost;
+ isux = isux + 1;
+ jsux = jsux + 1;
+ else
+ x2(irun,:) = ix2(b,:);
+ logpo2(irun) = ilogpo2(b);
+ end
+ prtfrc = j/nruns(b);
+ if exist('OCTAVE_VERSION')
+ if mod(j, 10) == 0,
+ printf('MH: Computing Metropolis-Hastings (chain %d/%d): %3.f%% done, acception rate: %3.f%%\r', b, nblck, 100 * prtfrc, 100 * isux / j);
+ end
+ if mod(j,50)==0 & whoiam,
+% keyboard;
+ waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) '), ' sprintf('accept. %3.f%%%%', 100 * isux/j)];
+ fMessageStatus(prtfrc,whoiam,waitbarString, '', options_.parallel(ThisMatlab), MasterName, DyMo)
+ end
+ else
+ if mod(j, 3)==0 & ~whoiam
+ waitbar(prtfrc,hh,[ '(' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)]);
+ elseif mod(j,50)==0 & whoiam,
+% keyboard;
+ waitbarString = [ '(' int2str(b) '/' int2str(options_.mh_nblck) ') ' sprintf('%f done, acceptation rate %f',prtfrc,isux/j)];
+ fMessageStatus(prtfrc,whoiam,waitbarString, waitbarTitle, options_.parallel(ThisMatlab), MasterName, DyMo)
+ end
+ end
+
+ if (irun == InitSizeArray(b)) | (j == nruns(b)) % Now I save the simulations
+ save([MhDirectoryName '/' ModelName '_mh' int2str(NewFile(b)) '_blck' int2str(b) '.mat'],'x2','logpo2');
+ fidlog = fopen([MhDirectoryName '/metropolis.log'],'a');
+ fprintf(fidlog,['\n']);
+ fprintf(fidlog,['%% Mh' int2str(NewFile(b)) 'Blck' int2str(b) ' (' datestr(now,0) ')\n']);
+ fprintf(fidlog,' \n');
+ fprintf(fidlog,[' Number of simulations.: ' int2str(length(logpo2)) '\n']);
+ fprintf(fidlog,[' Acceptation rate......: ' num2str(jsux/length(logpo2)) '\n']);
+ fprintf(fidlog,[' Posterior mean........:\n']);
+ for i=1:length(x2(1,:))
+ fprintf(fidlog,[' params:' int2str(i) ': ' num2str(mean(x2(:,i))) '\n']);
+ end
+ fprintf(fidlog,[' log2po:' num2str(mean(logpo2)) '\n']);
+ fprintf(fidlog,[' Minimum value.........:\n']);;
+ for i=1:length(x2(1,:))
+ fprintf(fidlog,[' params:' int2str(i) ': ' num2str(min(x2(:,i))) '\n']);
+ end
+ fprintf(fidlog,[' log2po:' num2str(min(logpo2)) '\n']);
+ fprintf(fidlog,[' Maximum value.........:\n']);
+ for i=1:length(x2(1,:))
+ fprintf(fidlog,[' params:' int2str(i) ': ' num2str(max(x2(:,i))) '\n']);
+ end
+ fprintf(fidlog,[' log2po:' num2str(max(logpo2)) '\n']);
+ fprintf(fidlog,' \n');
+ fclose(fidlog);
+ jsux = 0;
+ if j == nruns(b) % I record the last draw...
+ record.LastParameters(b,:) = x2(end,:);
+ record.LastLogLiK(b) = logpo2(end);
+ end
+ % size of next file in chain b
+ InitSizeArray(b) = min(nruns(b)-j,MAX_nruns);
+ % initialization of next file if necessary
+ if InitSizeArray(b)
+ x2 = zeros(InitSizeArray(b),npar);
+ logpo2 = zeros(InitSizeArray(b),1);
+ NewFile(b) = NewFile(b) + 1;
+ irun = 0;
+ end
+ end
+ j=j+1;
+ irun = irun + 1;
+ end% End of the simulations for one mh-block.
+ record.AcceptationRates(b) = isux/j;
+ if exist('OCTAVE_VERSION')
+ printf('\n');
+ diary on;
+ elseif ~whoiam
+ close(hh);
+ end
+ record.Seeds(b).Normal = randn('state');
+ record.Seeds(b).Unifor = rand('state');
+ OutputFileName(jloop,:) = {[MhDirectoryName,filesep], [ModelName '_mh*_blck' int2str(b) '.mat']};
+end% End of the loop over the mh-blocks.
+
+
+myoutput.record = record;
+myoutput.irun = irun;
+myoutput.NewFile = NewFile;
+myoutput.OutputFileName = OutputFileName;
+
+