* 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-bf33cf982152time-shift
stepan
parent
79125b8990
commit
1e4a7c2c71
|
|
@ -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
|
||||
|
|
@ -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;
|
||||
Loading…
Reference in New Issue