%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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_MC3_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\m8_MC3 betamat sig2mat gammat X y;
k=size(X,2);

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Model probabilities
%%%%%%%%%%%%%%%%%%%%%%%%%%
R=size(gammat,2); %#number of model draws

[uniquemat,c1,c2]=unique(gammat','rows'); %picks unique rows of gammat only; 
% this will be (number of unique models x k)
% only takes 10 sec or so

%uniquemat contains the unique models, in rows
%c1 shows the index of the first occurrence of each unique model
%c2 shows the index for each model, with repetitions, exactly in the order
%of occurrence in gammat

int=tabulate(c2);
uniquecount=int(:,2); %number of repetitions for unique models, same order as in uniquemat
U=length(uniquecount);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Rank models based on empirical frequency
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
int=[uniquemat uniquecount];
si=size(int,2);
int=sortrows(int,-si); %sorts all rows based on last column, 
% in descending order

Mmat=int(:,1:end-1);%already sorted
Mcount=int(:,end);%already sorted

relcount=Mcount/R; %relative frequency for each model

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

%Show UP TO top ten models
%%%%%%%%%%%%%%%%%%%%%%%%%%%
out=int(1:min(10,length(Mcount)),:)';
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');


relcount=Mcount/R; %relative frequency for each model

save c:\klaus\AAEC6564\mlab\worksp\m8_MC3_modProb Mmat Mcount relcount;


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




