%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2 component mixture model (2CMM), no covariates
% based on KPT section 15.5
% this version uses y generated by a mixture density
% here we will re-arrange the GS output sorted by magnitude of variance
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
tic;  % start stop watch

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

load c:\klaus\AAEC6564\mlab\worksp\2CMM_data; %to get the true values
y=ymix;
n=length(y);
G=2; %number of classes
k=size(X,2);

load c:\klaus\AAEC6564\mlab\worksp\mod10_2CMMmix beta1mat beta2mat ...
    sig11mat sig22mat cpmat Zmat;

% we want the larger variance draw to indicate regime 1
f1=find(sig11mat>sig22mat);
f2=find(sig11mat<sig22mat);

%make sure they are perfect complements
%tt=length(f1)+length(f2) % OK, 10,000, same as r2

beta1matNEW=beta1mat;
beta1matNEW(f2)=beta2mat(f2); %switch
beta2matNEW=beta2mat;
beta2matNEW(f2)=beta1mat(f2); %switch

sig11matNEW=sig11mat;
sig11matNEW(f2)=sig22mat(f2); %switch
sig22matNEW=sig22mat;
sig22matNEW(f2)=sig11mat(f2); %switch

cp1mat=cpmat(1,:);
cp2mat=cpmat(2,:);

cp1matNEW=cp1mat;
cp1matNEW(f2)=cp2mat(f2); %switch
cp2matNEW=cp2mat;
cp2matNEW(f2)=cp1mat(f2); %switch

cpmatNEW=[cp1matNEW;cp2matNEW];

'GS done'
% put all draws together & run diagnostics
allmat=[beta1matNEW;beta2matNEW;sig11matNEW;sig22matNEW;cpmatNEW];
kdiag=klausdiagnostics_greater0(allmat);

% Open log file 
[fid]=fopen('c:\Klaus\AAEC6564\mlab\logs\mod10_2CMMmix_rearranged.txt','w');
if fid==-1;
   warning('File could not be opened');
else;
   disp('File opened successfully');
end;

fprintf(fid,'true beta1\t%6.3f\n',beta1true'); 
fprintf(fid,'true beta2\t%6.3f\n',beta2true'); 
fprintf(fid,'\n');

sigAlltrue=[sig11true;sig22true];
fprintf(fid,'true sig2s \t%6.3f\n',sigAlltrue'); 
fprintf(fid,'\n');
fprintf(fid,'\n');

fprintf(fid,'true regime probs \t%6.3f\n',cptrue'); 
fprintf(fid,'\n');
fprintf(fid,'\n');

% beta stuff
out=kdiag(1:k,:)';
fprintf(fid,'Output table for beta1 \n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\t\tIEF\t\tm*\t\tcd\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');

out=kdiag(k+1:2*k,:)';
fprintf(fid,'Output table for beta2 \n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\t\tIEF\t\tm*\t\tcd\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');

out=kdiag(2*k+1:2*k+2,:)';
fprintf(fid,'Output table for sig2s \n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\t\tIEF\t\tm*\t\tcd\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');

out=kdiag(2*k+3:end,:)';
fprintf(fid,'Output table for cp \n\n');
fprintf(fid,'mean\t\tstd\t\tp(>0)\t\tnse\t\tIEF\t\tm*\t\tcd\n\n');
fprintf(fid,'%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n',out);
fprintf(fid,'\n');


finish = toc/60;
fprintf(fid,'Time elapsed in minutes \n\n');
fprintf(fid,'%6.3f\n',finish);

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