stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Removed check_presence_consecutive_MC_files in McMCDiagnostics. Cosmetic changes.

time-shift
Stéphane Adjemian (Scylla) 2013-11-22 11:28:24 +01:00
parent a65197f273
commit f14aed8f6a
2 changed files with 28 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

56
matlab/McMCDiagnostics.m
View File

@ -47,13 +47,30 @@ npar = npar + estim_params_.np ;
MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
load_last_mh_history_file(MetropolisFolder, ModelName);
NumberOfMcFilesPerBlock = record.LastFileNumber;
NumberOfMcFilesPerBlock = size(dir([MetropolisFolder ,filesep, ModelName '_mh*_blck1.mat']),1);
% Check that the declared number of blocks is consistent with informations saved in mh-history files.
if ~isequal(nblck,record.Nblck)
disp(['Estimation::mcmc::diagnostics: The number of MCMC chains you declared (' num2str(nblck) ') is inconsistent with the information available in the mh-history files (' num2str(record.Nblck) ' chains)!'])
disp([' I reset the number of MCMC chains to ' num2str(record.Nblck) '.'])
nblck = record.Nblck;
end
% check if all previous files are there for block 1
check_presence_consecutive_MC_files(MetropolisFolder,ModelName,1)
% check if all the mh files are available.
issue_an_error_message = 0;
for b = 1:nblck
nfiles = length(dir([MetropolisFolder ,filesep, ModelName '_mh*_blck' num2str(b) '.mat']));
if ~isequal(NumberOfMcFilesPerBlock,nfiles)
issue_an_error_message = 1;
disp(['Estimation::mcmc::diagnostics: The number of MCMC files in chain ' num2str(b) ' is ' num2str(nfiles) ' while the mh-history files indicate that we should have ' num2str(NumberOfMcFilesPerBlock) ' MCMC files per chain!'])
end
end
if issue_an_error_message
error('Estimation::mcmc::diagnostics: I cannot proceed because some MCMC files are missing. Check your MCMC files...')
end
if nblck == 1 % Brooks and Gelman tests need more than one block
if nblck == 1 % Brooks and Gelman tests need more than one block
convergence_diagnostics_geweke=zeros(npar,4+2*length(options_.convergence.geweke.taper_steps));
if any(options_.convergence.geweke.geweke_interval<0) || any(options_.convergence.geweke.geweke_interval>1) || length(any(options_.convergence.geweke.geweke_interval<0))~=2 ...
|| (options_.convergence.geweke.geweke_interval(2)-options_.convergence.geweke.geweke_interval(1)<0)
@ -104,17 +121,6 @@ if nblck == 1 % Brooks and Gelman tests need more than one block
return;
end
for blck = 2:nblck
tmp = size(dir([MetropolisFolder ,filesep, ModelName '_mh*_blck' int2str(blck) '.mat']),1);
if tmp~=NumberOfMcFilesPerBlock
disp(['McMCDiagnostics:: The number of mh files in chain ' int2str(blck) ' is ' int2str(tmp) ' while'])
disp([' the number of mh files in chain 1 is ' int2str(mcfiles) '!'])
error('The number of mh files has to be constant across chains!')
end
% check if all previous files are there for block blck
check_presence_consecutive_MC_files(MetropolisFolder,ModelName,blck)
end
PastDraws = sum(record.MhDraws,1);
LastFileNumber = PastDraws(2);
LastLineNumber = record.MhDraws(end,3);
@ -130,7 +136,7 @@ tmp = zeros(NumberOfDraws*nblck,3);
UDIAG = zeros(NumberOfLines,6,npar);
if NumberOfDraws < Origin
disp('MCMC Diagnostics :: The number of simulations is to small to compute the MCMC convergence diagnostics.')
disp('Estimation::mcmc::diagnostics: The number of simulations is too small to compute the MCMC convergence diagnostics.')
return
end
@ -141,7 +147,7 @@ if TeX
fprintf(fidTeX,' \n');
end
disp('MCMC Diagnostics: Univariate convergence diagnostic, Brooks and Gelman (1998):')
disp('Estimation::mcmc::diagnostics: Univariate convergence diagnostic, Brooks and Gelman (1998):')
% The mandatory variables for local/remote parallel
% computing are stored in localVars struct.
@ -324,7 +330,6 @@ MDIAG = zeros(NumberOfLines,6);
for b = 1:nblck
startline = 0;
for n = 1:NumberOfMcFilesPerBlock
%load([MetropolisFolder '/' mcfiles(n,1,b).name],'logpo2');
load([MetropolisFolder '/' ModelName '_mh',int2str(n),'_blck' int2str(b) '.mat'],'logpo2');
nlogpo2 = size(logpo2,1);
tmp((b-1)*NumberOfDraws+startline+(1:nlogpo2),1) = logpo2;
@ -409,17 +414,4 @@ if TeX
fprintf(fidTeX,'\n');
fprintf(fidTeX,'% End Of TeX file.');
fclose(fidTeX);
end
function check_presence_consecutive_MC_files(MetropolisFolder,fname,blck)
% check if all previous files are there
files=ls([MetropolisFolder ,filesep, fname '_mh*_blck' int2str(blck) '.mat']); %list files
right_string=files(:,length([fname '_mh'])+1:end); %cut off left part of filename
k = cell2mat(strfind(cellstr(right_string),['_blck' int2str(blck) '.mat'])); %find index of position after number
file_numbers=str2num(right_string(:,1:k-1)); %get file number
if ~isempty(file_numbers) && ...
sum(sort(file_numbers)-(min(file_numbers):max(file_numbers))')~=0
error(['There are MH draw files missing within chain ', int2str(blck)])
end
end

16
matlab/McMCDiagnostics_core.m
View File

@ -41,7 +41,7 @@ end
% Reshape 'myinputs' for local computation.
% In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by:
MhDirectoryName=myinputs.MhDirectoryName;
MetropolisFolder=myinputs.MetropolisFolder;%myinputs.MetropolisFolder;
nblck=myinputs.nblck;
NumberOfMcFilesPerBlock=myinputs.NumberOfMcFilesPerBlock;
Origin=myinputs.Origin;
@ -55,8 +55,8 @@ M_=myinputs.M_;
if whoiam
Parallel=myinputs.Parallel;
end
if ~exist('MhDirectoryName'),
MhDirectoryName = CheckPath('metropolis',M_.dname);
if ~exist('MetropolisFolder'),
MetropolisFolder = CheckPath('metropolis',M_.dname);
end
ALPHA = 0.2; % increase too much with the number of simulations.
@ -83,19 +83,11 @@ for j=fpar:npar,
for b = 1:nblck
startline = 0;
for n = 1:NumberOfMcFilesPerBlock
%load([MhDirectoryName '/' mcfiles(n,1,b).name],'x2');
load([MhDirectoryName '/' M_.fname '_mh',int2str(n),'_blck' int2str(b) ...
'.mat'],'x2');
load([MetropolisFolder '/' M_.fname '_mh',int2str(n),'_blck' int2str(b) '.mat'],'x2');
nx2 = size(x2,1);
tmp((b-1)*NumberOfDraws+startline+(1:nx2),1) = x2(:,j);
% clear x2;
startline = startline + nx2;
end
% $$$ %load([MhDirectoryName '/' mcfiles(NumberOfMcFilesPerBlock,1,b).name],'x2');
% $$$ load([MhDirectoryName '/' M_.fname '_mh',int2str(NumberOfMcFilesPerBlock),'_blck' int2str(b) '.mat'],'x2');
% $$$ tmp((b-1)*NumberOfDraws+startline+1:(b-1)*NumberOfDraws+MAX_nruns*(LastFileNumber-1)+LastLineNumber,1) = x2(:,j);
% $$$ clear x2;
% $$$ startline = startline + LastLineNumber;
end
tmp(:,2) = kron(transpose(1:nblck),ones(NumberOfDraws,1));
tmp(:,3) = kron(ones(nblck,1),time');