diff --git a/matlab/CheckPath.m b/matlab/CheckPath.m index 0e09700f5..41bd30263 100644 --- a/matlab/CheckPath.m +++ b/matlab/CheckPath.m @@ -4,6 +4,6 @@ global M_ DirectoryName = [ M_.dname '/' type ]; -if ~isdir(DirectoryName) - mkdir('.',DirectoryName); +if ~isdir(DirectoryName) + mkdir('.',DirectoryName); end \ No newline at end of file diff --git a/matlab/McMCDiagnostics.m b/matlab/McMCDiagnostics.m index 572e4f397..f5f8910f4 100644 --- a/matlab/McMCDiagnostics.m +++ b/matlab/McMCDiagnostics.m @@ -3,7 +3,8 @@ function McmcDiagnostic global options_ estim_params_ M_ -DirectoryName = CheckPath('Plots/Diagnostics'); +DirectoryName = CheckPath('Output'); +MhDirectoryName = CheckPath('metropolis'); TeX = options_.TeX; nblck = options_.mh_nblck; @@ -18,8 +19,15 @@ npar = npar + estim_params_.ncn; npar = npar + estim_params_.np ; MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8); -load(['./' M_.dname '/metropolis/' M_.fname '_mh_history.mat']) +load([MhDirectoryName '/' M_.fname '_mh_history.mat']) +mcfiles = []; +for blck = 1:nblck + mcfiles = cat(3,mcfiles,dir([MhDirectoryName '/' M_.fname '_mh*_blck' int2str(blck) '.mat'])); +end +NumberOfMcFilesPerBlock = size(mcfiles,1); + + PastDraws = sum(record.MhDraws,1); LastFileNumber = PastDraws(2); LastLineNumber = record.MhDraws(end,3); @@ -27,233 +35,51 @@ NumberOfDraws = PastDraws(1); Origin = 1000; StepSize = ceil((NumberOfDraws-Origin)/100);% So that the computational time does not -ALPHA = 0.2; % increase too much with the number of simulations. +ALPHA = 0.2; % increase too much with the number of simulations. time = 1:NumberOfDraws; xx = Origin:StepSize:NumberOfDraws; NumberOfLines = length(xx); 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.') return end + if TeX - fidTeX = fopen([M_.dname '\TeX\' 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):') for j=1:npar - fprintf(' Parameter %d... ',j); - for b = 1:nblck - startline = 0; - for n = 1:LastFileNumber-1 - eval(['load ' M_.dname '/metropolis/' M_.fname '_mh' int2str(n) '_blck' int2str(b)]); - clear logpo2 post2; - tmp((b-1)*NumberOfDraws+startline+1:(b-1)*NumberOfDraws+MAX_nruns*n,1) = x2(:,j); - clear x2; - startline = startline + MAX_nruns; - end - eval(['load ' M_.dname '/metropolis/' M_.fname '_mh' int2str(LastFileNumber) '_blck' int2str(b)]); - clear logpo2 post2; - 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) = temp(CSUP,1)-temp(CINF,1); - moyenne = mean(temp(:,1));%% Pooled mean. - UDIAG(ligne,3,j) = sum((temp(:,1)-moyenne).^2)/(nblck*n-1); - UDIAG(ligne,5,j) = sum(abs(temp(:,1)-moyenne).^3)/(nblck*n-1); - for i=1:nblck - pmet = temp(find(temp(:,2)==i)); - UDIAG(ligne,2,j) = UDIAG(ligne,2,j) + pmet(csup,1)-pmet(cinf,1); - moyenne = mean(pmet,1); %% Within mean. - UDIAG(ligne,4,j) = UDIAG(ligne,4,j) + sum((pmet(:,1)-moyenne).^2)/(n-1); - UDIAG(ligne,6,j) = UDIAG(ligne,6,j) + sum(abs(pmet(:,1)-moyenne).^3)/(n-1); - end - end - fprintf('Done! \n'); -end -UDIAG(:,[2 4 6],:) = UDIAG(:,[2 4 6],:)/nblck; -disp(' ') -clear pmet temp moyenne CSUP CINF csup cinf n linea iter tmp; -pages = floor(npar/3); -k = 0; -for i = 1:pages - h = figure('Name','MCMC univariate diagnostic (Brooks and Gelman,1998)'); - 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)]); - eval(['print -dpdf ' DirectoryName '/' M_.fname '_udiag' int2str(i)]); - saveas(h,[DirectoryName '/' M_.fname '_udiag' int2str(i) '.fig']); - 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 - h = figure('Name','MCMC univariate diagnostic (Brooks and Gelman, 1998)'); - 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)]); - eval(['print -dpdf ' DirectoryName '/' M_.fname '_udiag' int2str(pages+1)]); - saveas(h,[DirectoryName '/' M_.fname '_udiag' int2str(pages+1) '.fig']); - 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 -end % if reste > 0 -clear UDIAG; -%% -%% Multivariate diagnostic. -%% -if TeX - fidTeX = fopen([M_.dname '/TeX/' 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 + fprintf(' Parameter %d... ',j); + for b = 1:nblck startline = 0; - for n = 1:LastFileNumber-1 - eval(['load ' M_.dname '/metropolis/' M_.fname '_mh' int2str(n) '_blck' int2str(b)]); - clear x2 post2; - tmp((b-1)*NumberOfDraws+startline+1:(b-1)*NumberOfDraws+MAX_nruns*n,1) = logpo2; - startline = startline+MAX_nruns; + for n = 1:NumberOfMcFilesPerBlock-1 + %eval(['load ' M_.dname '/metropolis/' M_.fname '_mh' int2str(n) '_blck' int2str(b)]); + load([MhDirectoryName '/' mcfiles(n,1,b).name],'x2'); + %clear logpo2 post2; + tmp((b-1)*NumberOfDraws+startline+1:(b-1)*NumberOfDraws+MAX_nruns*n,1) = x2(:,j); + clear x2; + startline = startline + MAX_nruns; end - eval(['load ' M_.dname '/metropolis/' M_.fname '_mh' int2str(LastFileNumber) '_blck' int2str(b)]); - clear x2 post2; - tmp((b-1)*NumberOfDraws+startline+1:(b-1)*NumberOfDraws+ MAX_nruns*(LastFileNumber-1)+LastLineNumber,1) = logpo2; -end -clear logpo2; -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 + load([MhDirectoryName '/' mcfiles(NumberOfMcFilesPerBlock,1,b).name],'x2'); + % eval(['load ' M_.dname '/metropolis/' M_.fname '_mh' int2str(LastFileNumber) '_blck' int2str(b)]); + % clear logpo2 post2; + 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; @@ -262,65 +88,267 @@ for iter = Origin:StepSize:NumberOfDraws 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); + UDIAG(ligne,1,j) = 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); + UDIAG(ligne,3,j) = sum((temp(:,1)-moyenne).^2)/(nblck*n-1); + UDIAG(ligne,5,j) = 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); + pmet = temp(find(temp(:,2)==i)); + UDIAG(ligne,2,j) = UDIAG(ligne,2,j) + pmet(csup,1)-pmet(cinf,1); + moyenne = mean(pmet,1); %% Within mean. + UDIAG(ligne,4,j) = UDIAG(ligne,4,j) + sum((pmet(:,1)-moyenne).^2)/(n-1); + UDIAG(ligne,6,j) = UDIAG(ligne,6,j) + sum(abs(pmet(:,1)-moyenne).^3)/(n-1); end + end + fprintf('Done! \n'); +end +UDIAG(:,[2 4 6],:) = UDIAG(:,[2 4 6],:)/nblck; +disp(' ') +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)]); + eval(['print -dpdf ' DirectoryName '/' M_.fname '_udiag' int2str(i)]); + saveas(h,[DirectoryName '/' M_.fname '_udiag' int2str(i) '.fig']); + 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; + end + end + eval(['print -depsc2 ' DirectoryName '/' M_.fname '_udiag' int2str(pages+1)]); + eval(['print -dpdf ' DirectoryName '/' M_.fname '_udiag' int2str(pages+1)]); + saveas(h,[DirectoryName '/' M_.fname '_udiag' int2str(pages+1) '.fig']); + 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 +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 = []; +end +tmp = zeros(NumberOfDraws*nblck,3); +MDIAG = zeros(NumberOfLines,6); +for b = 1:nblck + startline = 0; + for n = 1:NumberOfMcFilesPerBlock-1 + % eval(['load ' M_.dname '/metropolis/' M_.fname '_mh' int2str(n) '_blck' int2str(b)]); + % clear x2 post2; + load([MhDirectoryName '/' mcfiles(n,1,b).name],'logpo2'); + tmp((b-1)*NumberOfDraws+startline+1:(b-1)*NumberOfDraws+MAX_nruns*n,1) = logpo2; + startline = startline+MAX_nruns; + end + load([MhDirectoryName '/' mcfiles(NumberOfMcFilesPerBlock,1,b).name],'logpo2'); + % eval(['load ' M_.dname '/metropolis/' M_.fname '_mh' int2str(LastFileNumber) '_blck' int2str(b)]); + % clear x2 post2; + tmp((b-1)*NumberOfDraws+startline+1:(b-1)*NumberOfDraws+ MAX_nruns*(LastFileNumber-1)+LastLineNumber,1) = logpo2; +end +clear logpo2; +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); + 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; -h = figure('Name','Multivatiate diagnostic'); +if options_.nograph + h = figure('Name','Multivatiate diagnostic','Visible','off'); +else + 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']); eval(['print -dpdf ' DirectoryName '/' M_.fname '_mdiag']); saveas(h,[DirectoryName '/' M_.fname '_mdiag.fig']); 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 diff --git a/matlab/metropolis.m b/matlab/metropolis.m index ac9dd31f0..67004f2d8 100644 --- a/matlab/metropolis.m +++ b/matlab/metropolis.m @@ -99,7 +99,7 @@ if options_.load_mh_file == 0 record.LastLogLiK = zeros(nblck,1); record.LastFileNumber = AnticipatedNumberOfFiles+1; record.LastLineNumber = AnticipatedNumberOfLinesInTheLastFile; - save([DirectoryName '/' M_.fname '_mh_history'],'record'); + save([MhDirectoryName '/' M_.fname '_mh_history'],'record'); elseif options_.load_mh_file == 1% Here we consider previous mh files (previous mh did not crash). disp('MH: I''m loading past metropolis-hastings simulations...') file = dir([ MhDirectoryName '/' M_.fname '_mh_history.mat' ]);