%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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;