%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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_growth_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_growth betamat sig2mat gammat X y;
k=size(X,2); %X does not include a constant term
kf = k; %total number of "flexible coefficients"
%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
int=[int relcount];%
fprintf(fid,'total model space =\t%6.0f \n',2^kf);
fprintf(fid,'number of visited models =\t%6.0f \n',U);
fprintf(fid,'fraction of possible models visited =\t%6.3f \n',U/(2^kf));
fprintf(fid,'\n');
%Show UP TO top ten models
%%%%%%%%%%%%%%%%%%%%%%%%%%%
out=[Mcount(1:100) relcount(1:100)]';
fprintf(fid,'Most probable model: counts and emp. frequencies\n\n');
fprintf(fid,'%6.0f\t%6.8f\n',out);
fprintf(fid,'\n');
relcount=Mcount/R; %relative frequency for each model
save c:\klaus\AAEC6564\mlab\worksp\m8_MC3_growth_modProb Mmat Mcount relcount;
st=fclose(fid);
if st==0;
disp('File closed successfully');
else;
warning('Problem with closing file');
end;