diff mc_global_info1.m @ 18:062d46712995 tip

Moved mc_global_info1 back to public folder

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mc_global_info1.m	Mon Apr 02 21:50:43 2012 +0100
@@ -0,0 +1,78 @@
+% 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_fixpt2(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
+
+function p=mc_fixpt2(Pt,tol)
+	[U,S,V]=svd(Pt-eye(size(Pt)),0);
+	p=V(:,diag(S)<=tol);
+	p=max(0,p*diag(1./sum(p)));
+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)
+	[V,D]=eig(Pt);
+	p=V(:,abs(diag(D)-1)<tol);
+	p=real(p*diag(1./sum(p)));
+	p(abs(p)<eps)=0;
+end