samer@0: function Info=mc_global_info(T,pz)
samer@0: % mc_global_info - Three global information measures about stationary MC
samer@0: %
samer@0: % mc_global_info :: [[N,N]] ~'transmat' -> [[3]].
samer@0: %
samer@0: % mc_global_info :: 
samer@0: %    [[N,N]] ~'transmat',
samer@0: %    [[N]] ~'stationary distribution if already known'
samer@0: % -> [[3]].
samer@0: %
samer@0: % The three measures in order are
samer@0: %    Entropy rate H(X|Z)
samer@0: %    Redundancy I(X,Z)
samer@0: %    Predictive information rate I(X,Y|Z)
samer@7: % All in NATS, not bits.
samer@0: 
samer@0: 
samer@0: n=size(T,1);
samer@0: if nargin<2, pz=mc_fixpt(T); end
samer@0: TlogT=T.*slog(T);
samer@0: HXZ=-sum(TlogT,1)*pz;
samer@0: 
samer@0: pxz=T; % p(x|z)
samer@0: pyxz=repmat(T,[1,1,n]).*repmat(shiftdim(T,-1),[n,1,1]); % p(y,x|z)
samer@0: pyz=sum(pyxz,2); % p(y|z)
samer@0: pyz(pyz==0)=realmin;
samer@0: pxyz=permute(pyxz./repmat(pyz,[1,n,1]),[2,1,3]); % p(x|y,z)
samer@0: 
samer@0: HXYz=-shiftdim(sum(sum(permute(pyxz,[2,1,3]).*slog(pxyz),1),2),2); % H(X|Y,z)
samer@0: IXYZ=(-HXYz' - sum(pxz.*slog(pxz)))*pz;
samer@0: 
samer@0: IXZ=max(0,(sum(TlogT)-slog(pz)')*pz); 
samer@0: 
samer@0: 
samer@0: Info=[HXZ;IXZ;IXYZ];
samer@0: if any(~isreal(Info))
samer@0: 	disp('unreal information');
samer@0: 	keyboard
samer@0: elseif any(Info<0)
samer@0: 	if any(Info<-eps)
samer@0: 		disp('negative information');
samer@0: 		disp(Info);
samer@0: 		keyboard
samer@0: 	else
samer@0: 		Info=max(0,Info);
samer@0: 	end
samer@0: end
samer@0: 
samer@0: function y=slog(x)
samer@0: % slog - safe log, x<=0 => slog(x)=0
samer@0: %
samer@0: % slog :: real -> real.
samer@0: 
samer@0: 	x(x<=0)=1;
samer@0: 	y=log(x);
samer@0: