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