%Eigenvalues for every 4 hours 
ind3=5; 
for j=1:645        
i1=ind3+(j-1);i2=ind3+29+(j-1);%jmax=645
%i1=ind3+15*(j-1);i2=ind3+29+15*(j-1);%jmax=43;  
%i1=ind3+18*(j-1);i2=ind3+35+18*(j-1);%jmax=36; 
%i1=ind3+72*(j-1);i2=ind3+76+72*(j-1);%jmax=8;
Mdh=Mdreq(:,i1:i2);R=cov(Mdh');[V,D]=eig(R);
ceigv=-(sort(-diag(D)));
sece(j)=ceigv(2)/ceigv(1);
%figure(1),plot(ceigv/max(ceigv)) 

%Sphericity Test 
p=26;N=i2-i1;
for k=0:p-1   
cdummy=ceigv(k+1:p).^( 1/(p-k) );
a=prod(cdummy);b=sum( ceigv(k+1:p) )/(p-k);
aic(j,k+1)=-2*(p-k)*N*log(a/b)+2*k*(2*p-k);
mdl(j,k+1)=-(p-k)*N*log(a/b)+0.5*k*(2*p-k)*log(N);
end 
end 
figure(2),mesh(aic),title('AIC estimator'),shading interp,view(2)
colormap('mastcolors'),colorbar 
figure(3),mesh(mdl),title('MDL estimator'),shading interp,view(2)  
colormap('mastcolors'),colorbar 
[dummy,raic]=min(aic');figure(4),plot(raic),title('Minimum index for AIC')
[dummy,rmdl]=min(mdl');figure(5),plot(rmdl),title('Minimum index for MDL')

%Sphericity Test for the all data 
%p=26;N=max(size(Mdreq));
%R=cov(Mdreq');[V D]=eig(R);ceigv=-sort(-diag(D));  
%for k=0:p-1
%cdummy=ceigv(k+1:p).^( 1/(p-k) );
%a=prod(cdummy);b=sum( ceigv(k+1:p) )/(p-k);
%taic(k+1)=-2*(p-k)*N*log(a/b)+2*k*(2*p-k);
%tmdl(k+1)=-(p-k)*N*log(a/b)+0.5*k*(2*p-k);%*log(N);
%end

