% Matlab code for demonstarting cross-validation with BUGS
% 
% Author: Aki Vehtari <Aki.Vehtari@hut.fi>

% Make indeces for cross-validation. Here we use n-fold-CV because
% it does not take much more time than 10-fold-CV, if n would be
% large, it would be trivial to change this to 10-fold-CV
ii=1:21;
for i=1:21
  itst{i}=ii(i);
  itr{i}=setdiff(ii,itst{i});
end

% This is used to simplify pairwise comparison of the models
wwc=dirrand(n,1000);

% load data
d=load('stacks.dat');
[n,p]=size(d);p=p-1;
% reserve memory for arrays
lps=zeros(1000,N);
lpss=zeros(1000,N,10);
% make cross-validation for all 5 models
for mi=1:5
  fprintf('%d: ',mi);
  % First make model with full data
  % (to get full posterior and to make bias correction)
  mat2bugs('stacks_dat.txt','Y',d(:,1),'x',d(:,2:end),'p',p,'N',n);
  unix(sprintf('echo "check(''Y:/cv_vs_dic/stacks/WinBUGS/stacks%d_bug.txt'')" > stacks_scr.txt ',mi));
  unix('cat stackstr_scr.txt >> stacks_scr.txt');
  unix(['wine /proj/ave/c/WinBUGS/WinBUGS14.exe /PAR "Y:/cv_vs_dic/stacks/WinBUGS/stacks_scr.txt" 2> /dev/null' ]);
  S=bugs2mat('bugsIndex.txt','bugs1.txt');
  ds=S.ds;     % BUGS computed predictive deviances
  trlps=ds/-2; % We like to handle predictive likelihoods
  for i1=1:length(itr)
    fprintf('.');
    % And CV models
    mat2bugs('stackscv_dat.txt','Y',d(itr{i1},1),'x',d(itr{i1},2:end),'Y2',d(itst{i1},1),1,'x2',d(itst{i1},2:end),'p',p,'N',length(itr{i1}),'N2',length(itst{i1}));
    unix(sprintf('echo "check(''Y:/cv_vs_dic/stacks/WinBUGS/stacks%dcv_bug.txt'')" > stacks_scr.txt ',mi));
    unix('cat stackscv_scr.txt >> stacks_scr.txt');
    unix(['wine /proj/ave/c/WinBUGS/WinBUGS14.exe /PAR "Y:/cv_vs_dic/stacks/WinBUGS/stacks_scr.txt" 2> /dev/null']);
    S=bugs2mat('bugsIndex.txt','bugs1.txt');
    lps(:,itst{i1})=S.ds2/-2;    % CV predictive
    lpss(:,itr{i1},i1)=S.ds/-2;  % For bias correction
    lpss(:,itst{i1},i1)=S.ds2/-2;% For bias correction
  end
  fprintf('\n');
  % Mean of the expected utility
  % for posterior predictive (train)
  trlC=sum(log(mean(exp(trlps))));
  m(mi).trlC=trlC;
  % for CV-predictive
  cvlC=sum(log(mean(exp(lps))));
  m(mi).cvlC=cvlC;
  % bias corrected
  cvclC=sum(log(mean(exp(lps))))+sum(log(mean(exp(trlps))))-mean(sum(log(mean(exp(lpss)))));
  m(mi).cvclC=cvclC;
  % corresponding DIC and peff
  m(mi).cvcD=-2*cvclC;
  m(mi).cvcp=trlC-cvclC;
  % reserve memory for arrays
  cvlCs=zeros(1000,1);
  cvclCs=zeros(1000,1);
  cvcps=zeros(1000,1);
  % precompute these outside of loop
  etrlps=exp(trlps);
  elps=exp(lps);
  elpss=exp(lpss);
  % sample from the distribution of the expected utility estimate
  % we use Bayesian Bootstrap (Rubin, 1981), that is,
  % we assume that the posterior probabilities for the samples 
  % of a random variable have a Dirichlet distribution
  for i2=1:1000
    % uncertainty from mcmc sampling
    w=dirrand(1000,n);
    lcpo=log(sum(elps.*w))';
    lcpotr=log(sum(etrlps.*w))';
    lcpocvtr=squeeze(log(sum(elpss.*repmat(w,[1 1 n]))));
    % uncertainty from not knowing the future data distribution
    % here we use wwc so that we have same random numbers for each model
    % this to makes pairwise comparison easier
    cvlCs(i2)=n*sum(lcpo.*wwc(:,i2));
    cvclCs(i2)=n*sum(lcpo.*wwc(:,i2))+n*sum(lcpotr.*wwc(:,i2))-mean(n*sum(lcpocvtr.*repmat(wwc(:,i2),[1 n])));
  end
  m(mi).cvlCs=cvlCs;
  m(mi).cvclCs=cvclCs;
end
fprintf('\n');
save m m

% Figures
lab={'1) normal','2) double exp.','3) logistic','4) {\it{}t}_4','5) {\it{}t}_4 as scale mixt.'};
set(0,'DefaultAxesLineWidth',1)
set(0,'DefaultLineLineWidth',1)

load ii
load m
n=21;

% Figure 1: CV vs DIC - predictive deviance
clf
a4
figscale=[1 1/2];
set(gcf,'pos',[27.4 3.8 10.795.*figscale(1) 15.26.*figscale(2)]);
set(gcf,'paperposition',[5.1 6.9 10.795.*figscale(1) 15.26.*figscale(2)]);
set(gcf,'DefaultAxesFontSize',13)
set(gcf,'DefaultTextFontSize',13)

clf
set(gca,'Box','on');
set(gca,'Position',[0.375 0.16 0.585 0.835])
ddic=[115.2 113.5 114.8 114.2 109.7];
dcvc=-2*[m(:).cvclC];
h(1)=line(dcvc,5:-1:1,'Color','b','Linewidth',3,'Marker','o','MarkerFaceColor','b','MarkerSize',10,'LineStyle','none');
h(2)=line(ddic,5:-1:1,'Color','r','Linewidth',3,'Marker','d','MarkerFaceColor','r','MarkerSize',10,'LineStyle','none');
set(gca,'YLim',[0.5 5.5],'YTick',[])
set(gca,'XGrid','on')
set(gca,'XLim',[106 132])
tx=min(get(gca,'XLim'))-diff(get(gca,'XLim'))/1.55;
for qi=1:5
  text(tx,6-qi,sprintf('%s',lab{qi}),'HorizontalAlignment','left')
end
h=legend([h],'CV','DIC');
xlabel('Expected predictive deviance')
print -depsc2 stack_dic.eps

% Figure 2: CV vs DIC - p_eff
clf
a4
figscale=[1 1/2];
set(gcf,'pos',[27.4 3.8 10.795.*figscale(1) 15.26.*figscale(2)]);
set(gcf,'paperposition',[5.1 6.9 10.795.*figscale(1) 15.26.*figscale(2)]);
set(gcf,'DefaultAxesFontSize',13)
set(gcf,'DefaultTextFontSize',13)

clf
set(gca,'Box','on');
set(gca,'Position',[0.335 0.14 0.635 0.85])
set(gca,'Position',[0.375 0.16 0.585 0.835])
pdic=[5.1 5.6 5.3 5.5 7.6];
pcvc=[m(:).cvcp];
h(1)=line(pcvc,5:-1:1,'Color','b','Linewidth',3,'Marker','o','MarkerFaceColor','b','MarkerSize',10,'LineStyle','none');
h(2)=line(pdic,5:-1:1,'Color','r','Linewidth',3,'Marker','d','MarkerFaceColor','r','MarkerSize',10,'LineStyle','none');
set(gca,'YLim',[0.5 5.5],'YTick',[])
set(gca,'XGrid','on')
set(gca,'XLim',[4 10])
tx=min(get(gca,'XLim'))-diff(get(gca,'XLim'))/1.55;
for qi=1:5
  text(tx,6-qi,sprintf('%s',lab{qi}),'HorizontalAlignment','left')
end
legend([h],'CV','DIC')
xlabel('Effective number of parameters')
print -depsc2 stack_peff.eps

% Figure 3: Pairwise comparison of models
% Compute pairwise comparison of t_4 scale mixture model to others
for i1=1:4
  q=-2*([m(5).cvclCs-m(i1).cvclCs]);
  qx=min(q)-0.5:0.02:max(q)+1;
  p=kernelp(q,qx);
  qxs{i1}=qx;
  qps{i1}=p;
end
clf
a4
figscale=[1 1/2];
set(gcf,'pos',[27.4 3.8 10.795.*figscale(1) 15.26.*figscale(2)]);
set(gcf,'paperposition',[5.1 6.9 10.795.*figscale(1) 15.26.*figscale(2)]);
set(gcf,'DefaultAxesFontSize',13)
set(gcf,'DefaultTextFontSize',13)
set(gca,'Box','on');
set(gca,'Position',[0.05 0.17 0.9 0.815])
qc=get(0,'DefaultAxesColorOrder');
for i1=1:4
  line(qxs{i1},qps{i1},'Color',qc(i1,:),'Linewidth',3)
end
set(gca,'YLim',[0 0.88],'Ytick',[],'XLim',[-10 10])
line([0 0],[0 max(get(gca,'YLim'))],'LineStyle','--','Color','k')
xlabel('Pairwise comparison of {\it{}t}_4 scale mixture to others')
h=legend(lab{1:4});
set(h,'pos', [0.51 0.6437 0.44 0.3182])
print -depsc2 stack_pcomp.eps