view private/mc_global_info.m @ 18:062d46712995 tip

Moved mc_global_info1 back to public folder

function [Info,ergodic]=mc_global_info(T,pz)
% mc_global_info - Three global information measures about stationary MC
%
% mc_global_info :: [[N,N]] ~'transmat' -> [[3]].
%
% mc_global_info :: 
%    [[N,N]] ~'transmat',
%    [[N]] ~'stationary distribution if already known'
% -> [[3]],bool.
%
% The three measures in order are
%    Entropy rate H(X|Z)
%    Redundancy I(X,Z)
%    Predictive information rate I(X,Y|Z)
% All in NATS, not bits.


n=size(T,1);
if nargin<2, pz=mc_fixpt(T,0.001); end

ergodic=(size(pz,2)==1);
if ergodic && all(isfinite(pz))

	TlogT=T.*slog(T);
	HXZ=-sum(TlogT,1)*pz;

	pxz=T; % p(x|z)
	pyxz=repmat(T,[1,1,n]).*repmat(shiftdim(T,-1),[n,1,1]); % p(y,x|z)
	pyz=sum(pyxz,2); % p(y|z)
	pyz(pyz==0)=realmin;
	pxyz=permute(pyxz./repmat(pyz,[1,n,1]),[2,1,3]); % p(x|y,z)

	HXYz=-shiftdim(sum(sum(permute(pyxz,[2,1,3]).*slog(pxyz),1),2),2); % H(X|Y,z)
	IXYZ=(-HXYz' - sum(pxz.*slog(pxz)))*pz;

	IXZ=max(0,(sum(TlogT)-slog(pz)')*pz); 


	Info=[HXZ;IXZ;IXYZ];
	if any(~isreal(Info))
		disp('unreal information');
		ergodic=0;
	elseif any(Info<0)
		if any(Info<-eps)
			disp('negative information');
			disp(pz');
			disp(Info');
			ergodic=0;
		else
			Info=max(0,Info);
		end
	end
else
	Info=[0;0;0];
end

function y=slog(x)
% slog - safe log, x<=0 => slog(x)=0
%
% slog :: real -> real.

	x(x<=0)=1;
	y=log(x);