view private/mc_global_info1.m @ 14:1d8ae1c1ee2f

test THEN commit...

% mc_global_info1 - Three global information measures about stationary MC
%
% mc_global_info1 :: 
%    [[N,N]] ~'transmat',
%    nonneg  ~'tolerance for fixed point check'
% -> [[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.
function [Info,ergodic]=mc_global_info1(T,tol)

	n=size(T,1);
	pz=mc_fixpt(T,tol); 
	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
end

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

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

% mc_fixpt - fixed point(s) of Markov chain (stationary state distribution)
%
% mc_fixpt :: [[K,K]]~'transition matrix'-> [[K,L]].
function p=mc_fixpt(Pt,tol)
	% use eigenvalues of transition matrix.
	[V,D]=eig(Pt);

	% find eigenvalues near 1 to within tolerance
	k=find(abs(diag(D)-1)<tol);
	
	Vk=V(:,k);
	%sumk=sum(V(:,k));
	%Vk=V(:,k)*diag(1./sumk); % normalise to unit sum
	%poss=abs(sumk)>0; % find those with non-zero sum
	%Vk=V(:,k(poss))*diag(1./sumk(poss)); % normalise to unit sum
	%p=Vk(:,all(Vk >= -eps)); % eigs with all positive components to within eps.

	%p=max(0,real(Vk*diag(1./sum(Vk))));
	p=real(Vk*diag(1./sum(Vk)));
	p(abs(p)<tol/10)=0;
end