From 1e4a7c2c717cf3818feb5ce79e11bb88db63f90c Mon Sep 17 00:00:00 2001 From: stepan Date: Thu, 24 Sep 2009 12:48:28 +0000 Subject: [PATCH] * Cosmetic changes. * Bug fix (with octave-3.2.2) in McMCDiagnostics.m git-svn-id: https://www.dynare.org/svn/dynare/trunk@2973 ac1d8469-bf42-47a9-8791-bf33cf982152 --- matlab/McMCDiagnostics.m | 458 +++++++++++++++++----------------- matlab/McMCDiagnostics_core.m | 105 ++++---- 2 files changed, 280 insertions(+), 283 deletions(-) diff --git a/matlab/McMCDiagnostics.m b/matlab/McMCDiagnostics.m index df798878c..02c504efe 100644 --- a/matlab/McMCDiagnostics.m +++ b/matlab/McMCDiagnostics.m @@ -37,7 +37,7 @@ TeX = options_.TeX; nblck = options_.mh_nblck; % Brooks and Gelman tests need more than one block if nblck == 1 - return; + return; end npar = estim_params_.nvx; npar = npar + estim_params_.nvn; @@ -48,13 +48,16 @@ MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8); load([MhDirectoryName '/' M_.fname '_mh_history.mat']) -mcfiles = []; -for blck = 1:nblck - mcfiles = cat(3,mcfiles,dir([MhDirectoryName ,filesep, M_.fname '_mh*_blck' int2str(blck) '.mat'])); +NumberOfMcFilesPerBlock = size(dir([MhDirectoryName ,filesep, M_.fname '_mh*_blck1.mat']),1); +for blck = 2:nblck + tmp = size(dir([MhDirectoryName ,filesep, M_.fname '_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 end -NumberOfMcFilesPerBlock = size(mcfiles,1); - PastDraws = sum(record.MhDraws,1); LastFileNumber = PastDraws(2); LastLineNumber = record.MhDraws(end,3); @@ -75,13 +78,12 @@ if NumberOfDraws < Origin end if TeX - fidTeX = fopen([DirectoryName '/' M_.fname '_UnivariateDiagnostics.TeX'],'w'); - fprintf(fidTeX,'%% TeX eps-loader file generated by McmcDiagnostics.m (Dynare).\n'); - fprintf(fidTeX,['%% ' datestr(now,0) '\n']); - fprintf(fidTeX,' \n'); + fidTeX = fopen([DirectoryName '/' M_.fname '_UnivariateDiagnostics.TeX'],'w'); + fprintf(fidTeX,'%% TeX eps-loader file generated by McmcDiagnostics.m (Dynare).\n'); + fprintf(fidTeX,['%% ' datestr(now,0) '\n']); + fprintf(fidTeX,' \n'); end - disp('MCMC Diagnostics: Univariate convergence diagnostic, Brooks and Gelman (1998):') localVars.MhDirectoryName = MhDirectoryName; @@ -98,16 +100,12 @@ if isnumeric(options_.parallel),% | isunix, % for the moment exclude unix platfo fout = McMCDiagnostics_core(localVars,1,npar,0); UDIAG = fout.UDIAG; clear fout - else - [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, 1, npar,{},'McMCDiagnostics_core', localVars); - UDIAG = fout(1).UDIAG; for j=2:totCPU, UDIAG = cat(3,UDIAG ,fout(j).UDIAG); end - end % for j=1:npar @@ -162,175 +160,175 @@ clear pmet temp moyenne CSUP CINF csup cinf n linea iter tmp; pages = floor(npar/3); k = 0; for i = 1:pages - if options_.nograph - h = figure('Name','MCMC univariate diagnostic (Brooks and Gelman,1998)','Visible','off'); - else - h = figure('Name','MCMC univariate diagnostic (Brooks and Gelman,1998)'); - end - boxplot = 1; - if TeX - NAMES = []; - TEXNAMES = []; - end - for j = 1:3 % Loop over parameters - k = k+1; - [nam,namtex] = get_the_name(k,TeX); - for crit = 1:3% Loop over criteria - if crit == 1 - plt1 = UDIAG(:,1,k); - plt2 = UDIAG(:,2,k); - namnam = [nam , ' (Interval)']; - elseif crit == 2 - plt1 = UDIAG(:,3,k); - plt2 = UDIAG(:,4,k); - namnam = [nam , ' (m2)']; - elseif crit == 3 - plt1 = UDIAG(:,5,k); - plt2 = UDIAG(:,6,k); - namnam = [nam , ' (m3)']; - end - if TeX - NAMES = strvcat(NAMES,deblank(namnam)); - TEXNAMES = strvcat(TEXNAMES,deblank(namtex)); - end - subplot(3,3,boxplot); - plot(xx,plt1,'-b'); % Pooled - hold on; - plot(xx,plt2,'-r'); % Within (mean) - hold off; - xlim([xx(1) xx(NumberOfLines)]) - title(namnam,'Interpreter','none') - boxplot = boxplot + 1; - end - end - eval(['print -depsc2 ' DirectoryName '/' M_.fname '_udiag' int2str(i) '.eps']); - if ~exist('OCTAVE_VERSION') - eval(['print -dpdf ' DirectoryName '/' M_.fname '_udiag' int2str(i)]); - end - if options_.nograph, set(h,'visible','on'), end - if ~exist('OCTAVE_VERSION') - saveas(h,[DirectoryName '/' M_.fname '_udiag' int2str(i) '.fig']); - end - if options_.nograph, close(h), end - if TeX - fprintf(fidTeX,'\\begin{figure}[H]\n'); - for jj = 1:size(NAMES,1) - fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TEXNAMES(jj,:))); + if options_.nograph + h = figure('Name','MCMC univariate diagnostic (Brooks and Gelman,1998)','Visible','off'); + else + h = figure('Name','MCMC univariate diagnostic (Brooks and Gelman,1998)'); end - fprintf(fidTeX,'\\centering \n'); - fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_udiag%s}\n',M_.fname,int2str(i)); - fprintf(fidTeX,'\\caption{Univariate convergence diagnostics for the Metropolis-Hastings.\n'); - fprintf(fidTeX,'The first, second and third columns are respectively the criteria based on\n'); - fprintf(fidTeX,'the eighty percent interval, the second and third moments.}'); - fprintf(fidTeX,'\\label{Fig:UnivariateDiagnostics:%s}\n',int2str(i)); - fprintf(fidTeX,'\\end{figure}\n'); - fprintf(fidTeX,'\n'); - end + boxplot = 1; + if TeX + NAMES = []; + TEXNAMES = []; + end + for j = 1:3 % Loop over parameters + k = k+1; + [nam,namtex] = get_the_name(k,TeX); + for crit = 1:3% Loop over criteria + if crit == 1 + plt1 = UDIAG(:,1,k); + plt2 = UDIAG(:,2,k); + namnam = [nam , ' (Interval)']; + elseif crit == 2 + plt1 = UDIAG(:,3,k); + plt2 = UDIAG(:,4,k); + namnam = [nam , ' (m2)']; + elseif crit == 3 + plt1 = UDIAG(:,5,k); + plt2 = UDIAG(:,6,k); + namnam = [nam , ' (m3)']; + end + if TeX + NAMES = strvcat(NAMES,deblank(namnam)); + TEXNAMES = strvcat(TEXNAMES,deblank(namtex)); + end + subplot(3,3,boxplot); + plot(xx,plt1,'-b'); % Pooled + hold on; + plot(xx,plt2,'-r'); % Within (mean) + hold off; + xlim([xx(1) xx(NumberOfLines)]) + title(namnam,'Interpreter','none') + boxplot = boxplot + 1; + end + end + eval(['print -depsc2 ' DirectoryName '/' M_.fname '_udiag' int2str(i) '.eps']); + if ~exist('OCTAVE_VERSION') + eval(['print -dpdf ' DirectoryName '/' M_.fname '_udiag' int2str(i)]); + end + if options_.nograph, set(h,'visible','on'), end + if ~exist('OCTAVE_VERSION') + saveas(h,[DirectoryName '/' M_.fname '_udiag' int2str(i) '.fig']); + end + if options_.nograph, close(h), end + if TeX + fprintf(fidTeX,'\\begin{figure}[H]\n'); + for jj = 1:size(NAMES,1) + fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TEXNAMES(jj,:))); + end + fprintf(fidTeX,'\\centering \n'); + fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_udiag%s}\n',M_.fname,int2str(i)); + fprintf(fidTeX,'\\caption{Univariate convergence diagnostics for the Metropolis-Hastings.\n'); + fprintf(fidTeX,'The first, second and third columns are respectively the criteria based on\n'); + fprintf(fidTeX,'the eighty percent interval, the second and third moments.}'); + fprintf(fidTeX,'\\label{Fig:UnivariateDiagnostics:%s}\n',int2str(i)); + fprintf(fidTeX,'\\end{figure}\n'); + fprintf(fidTeX,'\n'); + end end reste = npar-k; if reste - if reste == 1 - nr = 3; - nc = 1; - elseif reste == 2; - nr = 2; - nc = 3; - end - if TeX - NAMES = []; - TEXNAMES = []; - end - if options_.nograph - h = figure('Name','MCMC univariate diagnostic (Brooks and Gelman, 1998)','Visible','off'); - else - h = figure('Name','MCMC univariate diagnostic (Brooks and Gelman, 1998)'); - end - boxplot = 1; - for j = 1:reste - k = k+1; - [nam,namtex] = get_the_name(k,TeX); - for crit = 1:3 - if crit == 1 - plt1 = UDIAG(:,1,k); - plt2 = UDIAG(:,2,k); - namnam = [nam , ' (Interval)']; - elseif crit == 2 - plt1 = UDIAG(:,3,k); - plt2 = UDIAG(:,4,k); - namnam = [nam , ' (m2)']; - elseif crit == 3 - plt1 = UDIAG(:,5,k); - plt2 = UDIAG(:,6,k); - namnam = [nam , ' (m3)']; - end - if TeX - NAMES = strvcat(NAMES,deblank(namnam)); - TEXNAMES = strvcat(TEXNAMES,deblank(namtex)); - end - subplot(nr,nc,boxplot); - plot(xx,plt1,'-b'); % Pooled - hold on; - plot(xx,plt2,'-r'); % Within (mean) - hold off; - xlim([xx(1) xx(NumberOfLines)]); - title(namnam,'Interpreter','none'); - boxplot = boxplot + 1; + if reste == 1 + nr = 3; + nc = 1; + elseif reste == 2; + nr = 2; + nc = 3; end - end - eval(['print -depsc2 ' DirectoryName '/' M_.fname '_udiag' int2str(pages+1) '.eps']); - if ~exist('OCTAVE_VERSION') - eval(['print -dpdf ' DirectoryName '/' M_.fname '_udiag' int2str(pages+1)]); - end - if options_.nograph, set(h,'visible','on'), end - if ~exist('OCTAVE_VERSION') - saveas(h,[DirectoryName '/' M_.fname '_udiag' int2str(pages+1) '.fig']); - end - if options_.nograph, close(h), end - if TeX - fprintf(fidTeX,'\\begin{figure}[H]\n'); - for jj = 1:size(NAMES,1); - fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TEXNAMES(jj,:))); - end - fprintf(fidTeX,'\\centering \n'); - fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_udiag%s}\n',M_.fname,int2str(pages+1)); - if reste == 2 - fprintf(fidTeX,'\\caption{Univariate convergence diagnostics for the Metropolis-Hastings.\n'); - fprintf(fidTeX,'The first, second and third columns are respectively the criteria based on\n'); - fprintf(fidTeX,'the eighty percent interval, the second and third moments.}'); - elseif reste == 1 - fprintf(fidTeX,'\\caption{Univariate convergence diagnostics for the Metropolis-Hastings.\n'); - fprintf(fidTeX,'The first, second and third rows are respectively the criteria based on\n'); - fprintf(fidTeX,'the eighty percent interval, the second and third moments.}'); + if TeX + NAMES = []; + TEXNAMES = []; + end + if options_.nograph + h = figure('Name','MCMC univariate diagnostic (Brooks and Gelman, 1998)','Visible','off'); + else + h = figure('Name','MCMC univariate diagnostic (Brooks and Gelman, 1998)'); + end + boxplot = 1; + for j = 1:reste + k = k+1; + [nam,namtex] = get_the_name(k,TeX); + for crit = 1:3 + if crit == 1 + plt1 = UDIAG(:,1,k); + plt2 = UDIAG(:,2,k); + namnam = [nam , ' (Interval)']; + elseif crit == 2 + plt1 = UDIAG(:,3,k); + plt2 = UDIAG(:,4,k); + namnam = [nam , ' (m2)']; + elseif crit == 3 + plt1 = UDIAG(:,5,k); + plt2 = UDIAG(:,6,k); + namnam = [nam , ' (m3)']; + end + if TeX + NAMES = strvcat(NAMES,deblank(namnam)); + TEXNAMES = strvcat(TEXNAMES,deblank(namtex)); + end + subplot(nr,nc,boxplot); + plot(xx,plt1,'-b'); % Pooled + hold on; + plot(xx,plt2,'-r'); % Within (mean) + hold off; + xlim([xx(1) xx(NumberOfLines)]); + title(namnam,'Interpreter','none'); + boxplot = boxplot + 1; + end + end + eval(['print -depsc2 ' DirectoryName '/' M_.fname '_udiag' int2str(pages+1) '.eps']); + if ~exist('OCTAVE_VERSION') + eval(['print -dpdf ' DirectoryName '/' M_.fname '_udiag' int2str(pages+1)]); + end + if options_.nograph, set(h,'visible','on'), end + if ~exist('OCTAVE_VERSION') + saveas(h,[DirectoryName '/' M_.fname '_udiag' int2str(pages+1) '.fig']); + end + if options_.nograph, close(h), end + if TeX + fprintf(fidTeX,'\\begin{figure}[H]\n'); + for jj = 1:size(NAMES,1); + fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TEXNAMES(jj,:))); + end + fprintf(fidTeX,'\\centering \n'); + fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_udiag%s}\n',M_.fname,int2str(pages+1)); + if reste == 2 + fprintf(fidTeX,'\\caption{Univariate convergence diagnostics for the Metropolis-Hastings.\n'); + fprintf(fidTeX,'The first, second and third columns are respectively the criteria based on\n'); + fprintf(fidTeX,'the eighty percent interval, the second and third moments.}'); + elseif reste == 1 + fprintf(fidTeX,'\\caption{Univariate convergence diagnostics for the Metropolis-Hastings.\n'); + fprintf(fidTeX,'The first, second and third rows are respectively the criteria based on\n'); + fprintf(fidTeX,'the eighty percent interval, the second and third moments.}'); + end + fprintf(fidTeX,'\\label{Fig:UnivariateDiagnostics:%s}\n',int2str(pages+1)); + fprintf(fidTeX,'\\end{figure}\n'); + fprintf(fidTeX,'\n'); + fprintf(fidTeX,'% End Of TeX file.'); + fclose(fidTeX); end - fprintf(fidTeX,'\\label{Fig:UnivariateDiagnostics:%s}\n',int2str(pages+1)); - fprintf(fidTeX,'\\end{figure}\n'); - fprintf(fidTeX,'\n'); - fprintf(fidTeX,'% End Of TeX file.'); - fclose(fidTeX); - end end % if reste > 0 clear UDIAG; %% %% Multivariate diagnostic. %% if TeX - fidTeX = fopen([DirectoryName '/' M_.fname '_MultivariateDiagnostics.TeX'],'w'); - fprintf(fidTeX,'%% TeX eps-loader file generated by McmcDiagnostics.m (Dynare).\n'); - fprintf(fidTeX,['%% ' datestr(now,0) '\n']); - fprintf(fidTeX,' \n'); - NAMES = []; + fidTeX = fopen([DirectoryName '/' M_.fname '_MultivariateDiagnostics.TeX'],'w'); + fprintf(fidTeX,'%% TeX eps-loader file generated by McmcDiagnostics.m (Dynare).\n'); + fprintf(fidTeX,['%% ' datestr(now,0) '\n']); + fprintf(fidTeX,' \n'); + NAMES = []; end tmp = zeros(NumberOfDraws*nblck,3); MDIAG = zeros(NumberOfLines,6); for b = 1:nblck - startline = 0; - for n = 1:NumberOfMcFilesPerBlock - %load([MhDirectoryName '/' mcfiles(n,1,b).name],'logpo2'); - load([MhDirectoryName '/' M_.fname '_mh',int2str(n),'_blck' int2str(b) '.mat'],'logpo2'); - nlogpo2 = size(logpo2,1); - tmp((b-1)*NumberOfDraws+startline+(1:nlogpo2),1) = logpo2; - startline = startline+nlogpo2; - end + startline = 0; + for n = 1:NumberOfMcFilesPerBlock + %load([MhDirectoryName '/' mcfiles(n,1,b).name],'logpo2'); + load([MhDirectoryName '/' M_.fname '_mh',int2str(n),'_blck' int2str(b) '.mat'],'logpo2'); + nlogpo2 = size(logpo2,1); + tmp((b-1)*NumberOfDraws+startline+(1:nlogpo2),1) = logpo2; + startline = startline+nlogpo2; + end % $$$ %load([MhDirectoryName '/' mcfiles(NumberOfMcFilesPerBlock,1,b).name],'logpo2'); % $$$ load([MhDirectoryName '/' M_.fname '_mh',int2str(NumberOfMcFilesPerBlock),'_blck' int2str(b) '.mat'],'logpo2'); % $$$ tmp((b-1)*NumberOfDraws+startline+1:(b-1)*NumberOfDraws+ MAX_nruns*(LastFileNumber-1)+LastLineNumber,1) = logpo2; @@ -341,82 +339,82 @@ tmp(:,3) = kron(ones(nblck,1),time'); tmp = sortrows(tmp,1); ligne = 0; for iter = Origin:StepSize:NumberOfDraws - ligne = ligne+1; - linea = ceil(0.5*iter); - n = iter-linea+1; - cinf = round(n*ALPHA/2); - csup = round(n*(1-ALPHA/2)); - CINF = round(nblck*n*ALPHA/2); - CSUP = round(nblck*n*(1-ALPHA/2)); - temp = tmp(find((tmp(:,3)>=linea) & (tmp(:,3)<=iter)),1:2); - MDIAG(ligne,1) = temp(CSUP,1)-temp(CINF,1); - moyenne = mean(temp(:,1));%% Pooled mean. - MDIAG(ligne,3) = sum((temp(:,1)-moyenne).^2)/(nblck*n-1); - MDIAG(ligne,5) = sum(abs(temp(:,1)-moyenne).^3)/(nblck*n-1); - for i=1:nblck - pmet = temp(find(temp(:,2)==i)); - MDIAG(ligne,2) = MDIAG(ligne,2) + pmet(csup,1)-pmet(cinf,1); - moyenne = mean(pmet,1); %% Within mean. - MDIAG(ligne,4) = MDIAG(ligne,4) + sum((pmet(:,1)-moyenne).^2)/(n-1); - MDIAG(ligne,6) = MDIAG(ligne,6) + sum(abs(pmet(:,1)-moyenne).^3)/(n-1); - end + ligne = ligne+1; + linea = ceil(0.5*iter); + n = iter-linea+1; + cinf = round(n*ALPHA/2); + csup = round(n*(1-ALPHA/2)); + CINF = round(nblck*n*ALPHA/2); + CSUP = round(nblck*n*(1-ALPHA/2)); + temp = tmp(find((tmp(:,3)>=linea) & (tmp(:,3)<=iter)),1:2); + MDIAG(ligne,1) = temp(CSUP,1)-temp(CINF,1); + moyenne = mean(temp(:,1));%% Pooled mean. + MDIAG(ligne,3) = sum((temp(:,1)-moyenne).^2)/(nblck*n-1); + MDIAG(ligne,5) = sum(abs(temp(:,1)-moyenne).^3)/(nblck*n-1); + for i=1:nblck + pmet = temp(find(temp(:,2)==i)); + MDIAG(ligne,2) = MDIAG(ligne,2) + pmet(csup,1)-pmet(cinf,1); + moyenne = mean(pmet,1); %% Within mean. + MDIAG(ligne,4) = MDIAG(ligne,4) + sum((pmet(:,1)-moyenne).^2)/(n-1); + MDIAG(ligne,6) = MDIAG(ligne,6) + sum(abs(pmet(:,1)-moyenne).^3)/(n-1); + end end MDIAG(:,[2 4 6],:) = MDIAG(:,[2 4 6],:)/nblck; if options_.nograph - h = figure('Name','Multivatiate diagnostic','Visible','off'); + h = figure('Name','Multivatiate diagnostic','Visible','off'); else - h = figure('Name','Multivatiate diagnostic'); + h = figure('Name','Multivatiate diagnostic'); end boxplot = 1; for crit = 1:3 - if crit == 1 - plt1 = MDIAG(:,1); - plt2 = MDIAG(:,2); - namnam = 'Interval'; - elseif crit == 2 - plt1 = MDIAG(:,3); - plt2 = MDIAG(:,4); - namnam = 'm2'; - elseif crit == 3 - plt1 = MDIAG(:,5); - plt2 = MDIAG(:,6); - namnam = 'm3'; - end - if TeX - NAMES = strvcat(NAMES,namnam); - end - subplot(3,1,boxplot); - plot(xx,plt1,'-b'); % Pooled - hold on - plot(xx,plt2,'-r'); % Within (mean) - hold off - xlim([xx(1) xx(NumberOfLines)]) - title(namnam,'Interpreter','none'); - boxplot = boxplot + 1; + if crit == 1 + plt1 = MDIAG(:,1); + plt2 = MDIAG(:,2); + namnam = 'Interval'; + elseif crit == 2 + plt1 = MDIAG(:,3); + plt2 = MDIAG(:,4); + namnam = 'm2'; + elseif crit == 3 + plt1 = MDIAG(:,5); + plt2 = MDIAG(:,6); + namnam = 'm3'; + end + if TeX + NAMES = strvcat(NAMES,namnam); + end + subplot(3,1,boxplot); + plot(xx,plt1,'-b'); % Pooled + hold on + plot(xx,plt2,'-r'); % Within (mean) + hold off + xlim([xx(1) xx(NumberOfLines)]) + title(namnam,'Interpreter','none'); + boxplot = boxplot + 1; end eval(['print -depsc2 ' DirectoryName '/' M_.fname '_mdiag.eps']); if ~exist('OCTAVE_VERSION') - eval(['print -dpdf ' DirectoryName '/' M_.fname '_mdiag']); + eval(['print -dpdf ' DirectoryName '/' M_.fname '_mdiag']); end if options_.nograph, set(h,'visible','on'), end if ~exist('OCTAVE_VERSION') - saveas(h,[DirectoryName '/' M_.fname '_mdiag.fig']); + saveas(h,[DirectoryName '/' M_.fname '_mdiag.fig']); end if options_.nograph, close(h), end if TeX - fprintf(fidTeX,'\\begin{figure}[H]\n'); - for jj = 1:3 - fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),' '); - end - fprintf(fidTeX,'\\centering \n'); - fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_mdiag}\n',M_.fname); - fprintf(fidTeX,'\\caption{Multivariate convergence diagnostics for the Metropolis-Hastings.\n'); - fprintf(fidTeX,'The first, second and third rows are respectively the criteria based on\n'); - fprintf(fidTeX,'the eighty percent interval, the second and third moments. The different \n'); - fprintf(fidTeX,'parameters are aggregated using the posterior kernel.}'); - fprintf(fidTeX,'\\label{Fig:MultivariateDiagnostics}\n'); - fprintf(fidTeX,'\\end{figure}\n'); - fprintf(fidTeX,'\n'); - fprintf(fidTeX,'% End Of TeX file.'); - fclose(fidTeX); + fprintf(fidTeX,'\\begin{figure}[H]\n'); + for jj = 1:3 + fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),' '); + end + fprintf(fidTeX,'\\centering \n'); + fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_mdiag}\n',M_.fname); + fprintf(fidTeX,'\\caption{Multivariate convergence diagnostics for the Metropolis-Hastings.\n'); + fprintf(fidTeX,'The first, second and third rows are respectively the criteria based on\n'); + fprintf(fidTeX,'the eighty percent interval, the second and third moments. The different \n'); + fprintf(fidTeX,'parameters are aggregated using the posterior kernel.}'); + fprintf(fidTeX,'\\label{Fig:MultivariateDiagnostics}\n'); + fprintf(fidTeX,'\\end{figure}\n'); + fprintf(fidTeX,'\n'); + fprintf(fidTeX,'% End Of TeX file.'); + fclose(fidTeX); end \ No newline at end of file diff --git a/matlab/McMCDiagnostics_core.m b/matlab/McMCDiagnostics_core.m index 1dba7e643..0f1c7f01b 100644 --- a/matlab/McMCDiagnostics_core.m +++ b/matlab/McMCDiagnostics_core.m @@ -6,7 +6,7 @@ end struct2local(myinputs); if ~exist('MhDirectoryName'), - MhDirectoryName = CheckPath('metropolis'); + MhDirectoryName = CheckPath('metropolis'); end ALPHA = 0.2; % increase too much with the number of simulations. @@ -14,68 +14,67 @@ tmp = zeros(NumberOfDraws*nblck,3); UDIAG = zeros(NumberOfLines,6,npar-fpar+1); % keyboard; - if whoiam -% keyboard; - waitbarString = ['Please wait... McMCDiagnostics (' int2str(fpar) 'of' int2str(npar) ')...']; -% waitbarTitle=['McMCDiagnostics ',Parallel(ThisMatlab).PcName]; - if Parallel(ThisMatlab).Local, +if whoiam + % keyboard; + waitbarString = ['Please wait... McMCDiagnostics (' int2str(fpar) 'of' int2str(npar) ')...']; + % waitbarTitle=['McMCDiagnostics ',Parallel(ThisMatlab).PcName]; + if Parallel(ThisMatlab).Local, waitbarTitle=['Local ']; - else + else waitbarTitle=[Parallel(ThisMatlab).PcName]; - end - fMessageStatus(0,whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab), MasterName, DyMo); - - end + end + fMessageStatus(0,whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab), MasterName, DyMo); +end for j=fpar:npar, fprintf(' Parameter %d... ',j); - 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'); - nx2 = size(x2,1); - tmp((b-1)*NumberOfDraws+startline+(1:nx2),1) = x2(:,j); -% clear x2; - startline = startline + nx2; - end + 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'); + 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'); - tmp = sortrows(tmp,1); - ligne = 0; - for iter = Origin:StepSize:NumberOfDraws - ligne = ligne+1; - linea = ceil(0.5*iter); - n = iter-linea+1; - cinf = round(n*ALPHA/2); - csup = round(n*(1-ALPHA/2)); - CINF = round(nblck*n*ALPHA/2); - CSUP = round(nblck*n*(1-ALPHA/2)); - temp = tmp(find((tmp(:,3)>=linea) & (tmp(:,3)<=iter)),1:2); - UDIAG(ligne,1,j-fpar+1) = temp(CSUP,1)-temp(CINF,1); - moyenne = mean(temp(:,1));%% Pooled mean. - UDIAG(ligne,3,j-fpar+1) = sum((temp(:,1)-moyenne).^2)/(nblck*n-1); - UDIAG(ligne,5,j-fpar+1) = sum(abs(temp(:,1)-moyenne).^3)/(nblck*n-1); - for i=1:nblck - pmet = temp(find(temp(:,2)==i)); - UDIAG(ligne,2,j-fpar+1) = UDIAG(ligne,2,j-fpar+1) + pmet(csup,1)-pmet(cinf,1); - moyenne = mean(pmet,1); %% Within mean. - UDIAG(ligne,4,j-fpar+1) = UDIAG(ligne,4,j-fpar+1) + sum((pmet(:,1)-moyenne).^2)/(n-1); - UDIAG(ligne,6,j-fpar+1) = UDIAG(ligne,6,j-fpar+1) + sum(abs(pmet(:,1)-moyenne).^3)/(n-1); end - end - fprintf('Done! \n'); - if whoiam, -% keyboard; - waitbarString = [ 'Parameter ' int2str(j) '/' int2str(npar) ' done.']; - fMessageStatus((j-fpar+1)/(npar-fpar+1),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab), MasterName, DyMo) - end + tmp(:,2) = kron(transpose(1:nblck),ones(NumberOfDraws,1)); + tmp(:,3) = kron(ones(nblck,1),time'); + tmp = sortrows(tmp,1); + ligne = 0; + for iter = Origin:StepSize:NumberOfDraws + ligne = ligne+1; + linea = ceil(0.5*iter); + n = iter-linea+1; + cinf = round(n*ALPHA/2); + csup = round(n*(1-ALPHA/2)); + CINF = round(nblck*n*ALPHA/2); + CSUP = round(nblck*n*(1-ALPHA/2)); + temp = tmp(find((tmp(:,3)>=linea) & (tmp(:,3)<=iter)),1:2); + UDIAG(ligne,1,j-fpar+1) = temp(CSUP,1)-temp(CINF,1); + moyenne = mean(temp(:,1));%% Pooled mean. + UDIAG(ligne,3,j-fpar+1) = sum((temp(:,1)-moyenne).^2)/(nblck*n-1); + UDIAG(ligne,5,j-fpar+1) = sum(abs(temp(:,1)-moyenne).^3)/(nblck*n-1); + for i=1:nblck + pmet = temp(find(temp(:,2)==i)); + UDIAG(ligne,2,j-fpar+1) = UDIAG(ligne,2,j-fpar+1) + pmet(csup,1)-pmet(cinf,1); + moyenne = mean(pmet,1); %% Within mean. + UDIAG(ligne,4,j-fpar+1) = UDIAG(ligne,4,j-fpar+1) + sum((pmet(:,1)-moyenne).^2)/(n-1); + UDIAG(ligne,6,j-fpar+1) = UDIAG(ligne,6,j-fpar+1) + sum(abs(pmet(:,1)-moyenne).^3)/(n-1); + end + end + fprintf('Done! \n'); + if whoiam, + % keyboard; + waitbarString = [ 'Parameter ' int2str(j) '/' int2str(npar) ' done.']; + fMessageStatus((j-fpar+1)/(npar-fpar+1),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab), MasterName, DyMo) + end end myoutput.UDIAG = UDIAG; \ No newline at end of file