
% 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

  % note that the SSS models have x and y variables all standardised, and so
  % in the linear model (only) there is no intercept term. 

% -----------
% 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):(2+2p)       - posterior mode of the regression parameter vector beta 
    % elements (3+2p):(2+2p+p*p)  - estimated posterior variance matrix of beta in vectorised form 
    % final two elements:         - (s,d) - the residual SD estimate, and the posterior degees of freedom
%

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); 
my=mean(y(itrain)); sy=std(y(itrain)); Y=(y-my)/sy;

Fit=zeros(n,k);   sig=zeros(1,k); df=zeros(1,k);    % to save fits and predictions 
pmp=zeros(N,1);                                     % to save marginal inclusive probs
figure(1); clf
tau=1; del=3; 
aa=[min(y) max(y)];

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,2+p+(1:p))';                    % post mean of regn parameters
         A=X(ig,:)';                    	    % design matrix 
      else
          A=0;b=0;
      end
      sig(j)=models(j,end-1); df(j)=models(j,end);
      Fit(:,j)=my+sy*(A*b);                         % fitted & predicted linear regn
     % predictions for hold-out validation cases are already in the above
      clf; ii=ones(1,n); ii(ivalid)=0; 
      if (p>0)
          pairstv(Fit(:,j),y,itrain,ivalid)
          xlabel('fit');ylabel('data'); hold on; plot(aa,aa,'k:'); hold off
          title(['Model ',int2str(j),' : p=',int2str(p),' : Prob=',num2str(pm(j),2)])
      else
          show(y',ii,1)
      end
      fprintf('Hit a key to continue to next model ....\n') 
      pause
end

aveFit=Fit*pmk; clf
pairstv(aveFit,y,itrain,ivalid); title(['Model average'])
xlabel('aveFit');ylabel('data'); hold on; plot(aa,aa,'k:'); hold off
      
