samer@11: % mc_global_info1 - Three global information measures about stationary MC
samer@11: %
samer@11: % mc_global_info1 :: 
samer@11: %    [[N,N]] ~'transmat',
samer@11: %    nonneg  ~'tolerance for fixed point check'
samer@11: % -> [[3]],bool.
samer@11: %
samer@11: % The three measures in order are
samer@11: %    Entropy rate H(X|Z)
samer@11: %    Redundancy I(X,Z)
samer@11: %    Predictive information rate I(X,Y|Z)
samer@11: % All in NATS, not bits.
samer@14: function [Info,ergodic]=mc_global_info1(T,tol)
samer@11: 
samer@14: 	n=size(T,1);
samer@15: 	pz=mc_fixpt2(T,tol); 
samer@14: 	ergodic=(size(pz,2)==1);
samer@14: 	if ergodic && all(isfinite(pz))
samer@11: 
samer@14: 		TlogT=T.*slog(T);
samer@14: 		HXZ=-sum(TlogT,1)*pz;
samer@11: 
samer@14: 		pxz=T; % p(x|z)
samer@14: 		pyxz=repmat(T,[1,1,n]).*repmat(shiftdim(T,-1),[n,1,1]); % p(y,x|z)
samer@14: 		pyz=sum(pyxz,2); % p(y|z)
samer@14: 		pyz(pyz==0)=realmin;
samer@14: 		pxyz=permute(pyxz./repmat(pyz,[1,n,1]),[2,1,3]); % p(x|y,z)
samer@11: 
samer@14: 		HXYz=-shiftdim(sum(sum(permute(pyxz,[2,1,3]).*slog(pxyz),1),2),2); % H(X|Y,z)
samer@14: 		IXYZ=(-HXYz' - sum(pxz.*slog(pxz)))*pz;
samer@11: 
samer@14: 		IXZ=max(0,(sum(TlogT)-slog(pz)')*pz); 
samer@11: 
samer@11: 
samer@14: 		Info=[HXZ;IXZ;IXYZ];
samer@14: 		if any(~isreal(Info))
samer@14: 			disp('unreal information');
samer@11: 			ergodic=0;
samer@14: 		elseif any(Info<0)
samer@14: 			if any(Info<-eps)
samer@14: 				disp('negative information');
samer@14: 				disp(pz');
samer@14: 				disp(Info');
samer@14: 				ergodic=0;
samer@14: 			else
samer@14: 				Info=max(0,Info);
samer@14: 			end
samer@11: 		end
samer@14: 	else
samer@14: 		Info=[0;0;0];
samer@11: 	end
samer@11: end
samer@11: 
samer@11: function y=slog(x)
samer@11: % slog - safe log, x<=0 => slog(x)=0
samer@11: %
samer@11: % slog :: real -> real.
samer@11: 
samer@11: 	x(x<=0)=1;
samer@11: 	y=log(x);
samer@14: end
samer@11: 
samer@15: function p=mc_fixpt2(Pt,tol)
samer@15: 	[U,S,V]=svd(Pt-eye(size(Pt)),0);
samer@15: 	p=V(:,diag(S)<=tol);
samer@15: 	p=max(0,p*diag(1./sum(p)));
samer@15: end
samer@15: 
samer@15: 
samer@11: % mc_fixpt - fixed point(s) of Markov chain (stationary state distribution)
samer@11: %
samer@11: % mc_fixpt :: [[K,K]]~'transition matrix'-> [[K,L]].
samer@11: function p=mc_fixpt(Pt,tol)
samer@11: 	[V,D]=eig(Pt);
samer@15: 	p=V(:,abs(diag(D)-1)<tol);
samer@15: 	p=real(p*diag(1./sum(p)));
samer@15: 	p(abs(p)<eps)=0;
samer@11: end