%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get model probabilities and model-averaged results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rand('state',37); % set arbitrary seed for uniform draws
randn('state',37); % set arbitrary seed for normal draws

tic;  % start stop watch

[fid]=fopen('c:\klaus\AAEC6564\mlab\logs\mod8_SSVS_modProb.txt','w');
if fid==-1;
   warning('File could not be opened');
  return
else;
   disp('File opened successfully');
end;

load c:\klaus\AAEC6564\mlab\worksp\mod8_SSVS betamat sig2mat gammat X y;

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Model probabilities
%%%%%%%%%%%%%%%%%%%%%%%%%%
k=size(X,2);
Mmat = graycode(k-1); %lists all possible models as rows of zeros and 1
M=size(Mmat,1); %number of all possible models

% Get counter for visited models
% This works if your model space is not "huge" (a few 1000 is OK)
R=size(gammat,2); %#number of model draws
Mcount=zeros(M,1); %will collect empirical counts for visited models

for i=1:R
    int=repmat(gammat(:,i)',M,1); %J by 1
    int=abs(Mmat-int);
    f=find(sum(int,2)==0);
    Mcount(f)=Mcount(f)+1;
end
relcount=Mcount/R; 
%these are the empirical posterior model probabilities, out
% of all possible models
%[Mmat Mcount] OK, looks good
f=find(Mcount>0);
Mv=length(f);

fprintf(fid,'total model space =\t%6.0f \n',M);
fprintf(fid,'number of visited models =\t%6.0f \n',Mv);
fprintf(fid,'fraction of possible models visited =\t%6.3f \n',Mv/M);
fprintf(fid,'\n');


%Show UP TO top ten models
%%%%%%%%%%%%%%%%%%%%%%%%%%%
int=[Mmat relcount];
ik=size(int,2);
int=sortrows(int,-ik); %sort in descending order by last column (model count)
l=length(find(int(:,ik)>0));
out=int(1:min(10,l),:)';

fprintf(fid,'Most probable models \n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');


st=fclose(fid);
if st==0;
   disp('File closed successfully');
else;
   warning('Problem with closing file');
end;




