
% Look at variables selected and estimated coeffs of them in top k models, with
% predictions from each and the model average over top k models.  Also
% compute the corresponding conditional variable inclusion probabilities
% -- the latter in the vector pmp 

% -----------
% models array has following elements:  each row is one model of the 1:nm "top models"
% ordered in decreasing order of posterior probability. The number of columns is defined 
% by the largest model, and entries are NA/NaN for models smaller than the largest
% In each model/row the entries are: 
    % element  1:                 - dimension of the model = mumber of predictors p for this model
    % element  2:                 - log posterior probability of this model (the "score") 
    % elements 3:(2+p)            - the indices of the p variables in this model
    % elements (3+p):(3+2p)       - posterior mean of the regression parameter vector beta (includes intercept)
    % elements (4+2p):(4+4p+p*p)  - posterior variance matrix of beta (includes intercept) in vectorised form 
%

pmk=pm(1:k)/sum(pm(1:k));                           % condition on only the top k models chosen
X=(x-repmat(mean(x(:,itrain),2),1,n)) ./repmat(std(x(:,itrain)')',1,n); % needs since SSS standardises in model search

Fit=zeros(n,k); pFit=zeros(n,k);                    % to save fits and predictions 
pmp=zeros(N,1);                                     % to save marginal inclusive probs
figure(1); clf

for j=1:k
      p=models(j,1);                                % dim of this model
      if (p>0) 
          ig=models(j,2+(1:p));                     % predictors in model    
          pmp(ig)=pmp(ig)+pmk(j);         	        % posterior probs on predictors 
          b=models(j,3+p+(0:p))';                   % post mean of regn parameters
          A=[ones(n,1) X(ig,:)'];          	        % design matrix
      else
          b=models(j,3); A=[ones(n,1)]; 
      end    
      Fit(:,j)=A*b;                                 % fitted & predicted linear regn
      pFit(:,j)=1./(1+exp(-Fit(:,j)));              %    ... and plug-in probabilities
     % predictions for hold-out validation cases are already in the above
     showtv(pFit',y,itrain,ivalid,j); hold on; plot([0 n+1],[1 1]*sum(y)/n,'k:'); hold off
     title(['Model ',int2str(j),' : p=',int2str(p),' : Prob=',num2str(pm(j),2)])
	xlabel('Samples'); ylabel('Probabilities')
     if (length(ivalid)>0),
         legend('obs 0','obs 1','pred 0','pred 1','base','location','NorthEastOutside') 
     else
         legend('obs 0','obs 1','base','location','NorthEastOutside')
     end
     legend boxoff
     fprintf('Hit a key to continue to next model ....\n') 
     pause
end

avepFit=pFit*pmk;
clf; showtv(avepFit',y,itrain,ivalid,1); hold on; plot([0 n+1],[1 1]*sum(y)/n,'k:'); hold off
	xlabel('Samples'); ylabel('Probabilities')
	title(['Model average'])
     if (length(ivalid)>0),
         legend('obs 0','obs 1','pred 0','pred 1','base','location','NorthEastOutside') 
     else
         legend('obs 0','obs 1','base','location','NorthEastOutside')
     end
     legend boxoff
    
clear b A ig p a i 

