changeset 11:0e0f2805ef9c

Added new mechanism for checking Markov chains for uniqueness of stationary distribution;
new supporting files and new parameters to mt_init (see CHANGES).
Some functions require greater Matlab library, not included.

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/CHANGES	Sun Feb 26 23:11:10 2012 +0000
@@ -0,0 +1,17 @@
+Changed transition matrix sampling procedure to try to avoid Markov chains
+with mupltiple stationary distributions (fixed points).
+To this end:
+* mc_fixpt has a tolerance parameter - all eigenvalues within this much of
+  1 are considered fixed points.
+* mc_global_info1 (replaces mc_global_info) now checks that transition matrix
+  has only one fixed point and returns a flag if this is the case
+* sampling logic is in mt_resample - transition matrices not passing uniqueness
+  test are either resampled (method 1) or perturbed (method 2)
+* mt_init has new parameters, ErgMethod to select conditioning method and
+  tolerance to pass to mc_fixpt.
+
+
+Other changes 
+* Information measures are now displayed in bits, not nats.
+* Clicking on the title bar of a displayed transition matrix now displays
+  information about its eigenvalues and eigenvectors (mt_show_transmat.m).

--- a/mt_init.m	Wed Feb 22 21:24:44 2012 +0000
+++ b/mt_init.m	Sun Feb 26 23:11:10 2012 +0000
@@ -18,12 +18,19 @@
 %
 % The figure for scatter plots is fixed to figure 50 for now.
 
-function Sys=mt_init(A0,B0,A1,B1,L,Shuffle)
-	Sys.sample_transmats = @(k)sample_transmat_hdp(A0,B0,A1,B1,k,L);
+function Sys=mt_init(A0,B0,A1,B1,L,Shuffle,ErgMethod,Tol)
+	if nargin<7, ErgMethod=1; end
+	if nargin<6, Shuffle=0; end
+	if nargin<8, Tol=0.001; end
+	Sys.sample_transmats = @(k,l)sample_transmat_hdp(A0,B0,A1,B1,k,l);
 	Sys.transmats = {};
 	Sys.info      = {};
 	Sys.refpoints = [0,0;1,0;0,1]';
 	Sys.fig       = 50;
 	Sys.shuffle   = Shuffle;
+	Sys.ergmeth   = ErgMethod;
+	Sys.L         = L;
+	Sys.tol       = Tol;
 	Sys = mt_calibrate(Sys);
 end
+

--- a/mt_resample.m	Wed Feb 22 21:24:44 2012 +0000
+++ b/mt_resample.m	Sun Feb 26 23:11:10 2012 +0000
@@ -10,11 +10,42 @@
 % mt_init.
 
 function Sys=mt_resample(Sys,K)
-	[TT, II] = Sys.sample_transmats(K);
-	figure(Sys.fig); % 'name','info space'); 
+	figure(Sys.fig); 
+	cla;
+
+
+	L=Sys.L;
+	tol=Sys.tol;
+	II=zeros(L,3);
+	ok=zeros(1,L);
+
+	title(sprintf('sampling %d transition matrices...',L)); drawnow;
+	TT=Sys.sample_transmats(K,L);
+	title(sprintf('checking (tolerance=%f)...',tol)); drawnow;
+	for i=1:L, [II(i,:),ok(i)]=mc_global_info(TT(:,:,i)); end
+	todo=find(~ok);
+
+	if Sys.ergmeth==1
+		while ~isempty(todo)
+			title(sprintf('resampling %d nonergodic transition matrices...',length(todo)));
+			drawnow;
+			TT(:,:,todo)=Sys.sample_transmats(K,length(todo));
+			for i=todo, [II(i,:),ok(i)]=mc_global_info1(TT(:,:,i),tol); end
+			todo=find(~ok);
+		end
+	else
+		while ~isempty(todo)
+			title(sprintf('perturbing %d nonergodic transition matrices...',length(todo)));
+			drawnow;
+			for i=todo, j=randnat(K); TT(:,j,i)=stoch(TT(:,j,i)+0.1); end
+			for i=todo, [II(i,:),ok(i)]=mc_global_info1(TT(:,:,i),tol); end
+			todo=find(~ok);
+		end
+	end
+
 	Sys.transmats{K} = TT;
 	Sys.info{K}      = II;
-	Sys.hScat{K}     = scatc(II,II(:,3),16); 
+	Sys.hScat{K}     = scatc(II/log(2),II(:,3)/log(2),16); 
 	axis on; box on; grid off; 
 	lc=[0.4,0.4,0.4];
 	set(gca,'Color','none');
@@ -26,3 +57,7 @@
 	zlabel('pred-info rate');
 	rotate3d on;
 end
+
+
+
+

--- a/mt_show_transmat.m	Wed Feb 22 21:24:44 2012 +0000
+++ b/mt_show_transmat.m	Sun Feb 26 23:11:10 2012 +0000
@@ -4,11 +4,27 @@
 %    natural ~'voice id'
 %    [[K,K]] ~'transition matrix'
 % -> action void.
+%
+% Also displays global information measures in bits in the title.
 
 function mt_show_transmat(Id,T,I)
 	figure(Id); 
 	imagesc(T,[0,1]); 
 	axis xy;
 	if nargin<3, I=mc_global_info(T); end
-	title(sprintf('H=%.2f, R=%.2f, PI=%.2f',I(1),I(2),I(3)));
+	I=I/log(2); % convert to bits
+	htit=title(sprintf('H=%.2f, R=%.2f, PI=%.2f',I(1),I(2),I(3)));
+	addkbcallback(gcf);
+	addkbcallback(gcf,@(k)show_eigs(Id,T));
+	set(htit,'ButtonDownFcn',@(a,b)show_eigs(Id,T));
 end
+
+function show_eigs(Id,T) 
+	figure(40);
+	[V,D]=eig(T); d=diag(D);
+	k=find(abs(d-1)<0.001);
+	subplot(3,1,1); plot_cvec(diag(D)); 
+	title(sprintf('%d: eigenvalues %s',Id,mat2str(d(k))));
+	subplot(3,1,2); plot_cvec(sum(V)); title('sum eigenvectors');%ylim([-1,1]);
+	subplot(3,1,3); plotseq(@plot_cvec,window(V)); title('eigenvector');
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/private/addkbcallback.m	Sun Feb 26 23:11:10 2012 +0000
@@ -0,0 +1,19 @@
+function addkbcallback(fig,handler)
+% addkbcallback - Add keypress handler to figure or remove all handlers
+
+	if nargin>1,
+		cb=get(fig,'KeypressFcn');
+		if isempty(cb)
+			set(fig,'KeypressFcn',{@kpcb,{handler}});
+		else
+			set(fig,'KeypressFcn',{cb{1},[cb{2},{handler}]});
+		end
+	else
+		set(fig,'KeypressFcn',[]);
+	end
+end
+
+function kpcb(a,b,handlers), 
+	for i=1:length(handlers), h=handlers{i}; h(b); end; 
+end
+

--- a/private/mc_fixpt.m	Wed Feb 22 21:24:44 2012 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,34 +0,0 @@
-function p=mc_fixpt(Pt)
-% mc_fixpt - fixed point of Markov chain (stationary state distribution)
-%
-% mc_fixpt :: [[K,K]]~'transition matrix'-> [[K]].
-
-
-% for ergodic HMM we want to get stationary distribution 
-% use eigenvalues of transition matrix.
-[V,D]=eig(Pt);
-
-% one of the eigenvalues should be 1
-[dummy,k]=max(-abs(diag(D)-1));
-
-% check that it is really close
-if abs(D(k,k)-1)>0.001
-	disp('mc_fixpt: failed to find eigenvalue of 1 in ');
-	disp(Pt)
-	error('mc_fixpt:noeig','failed to find eigenvalue of 1');
-end
-
-% get eigenvector and flip it over if all negative
-v=V(:,k); if sum(v)<0, v=-v; end
-
-if ~isreal(v)
-	disp('mc_fixpt: discarding complex parts of eigenvector'); 
-	v=real(v);
-end
-if any(v<0)
-	%fprintf('mc_fixpt: discarding negative parts of eigenvector: %s\n',mat2str(v(v<0))); 
-	v=max(0,v);
-end
-
-% final normalisation
-p=v/sum(v);

--- a/private/mc_global_info.m	Wed Feb 22 21:24:44 2012 +0000
+++ b/private/mc_global_info.m	Sun Feb 26 23:11:10 2012 +0000
@@ -1,4 +1,4 @@
-function Info=mc_global_info(T,pz)
+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]].
@@ -6,7 +6,7 @@
 % mc_global_info :: 
 %    [[N,N]] ~'transmat',
 %    [[N]] ~'stationary distribution if already known'
-% -> [[3]].
+% -> [[3]],bool.
 %
 % The three measures in order are
 %    Entropy rate H(X|Z)
@@ -16,34 +16,42 @@
 
 
 n=size(T,1);
-if nargin<2, pz=mc_fixpt(T); end
-TlogT=T.*slog(T);
-HXZ=-sum(TlogT,1)*pz;
+if nargin<2, pz=mc_fixpt(T,0.001); end
 
-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)
+ergodic=(size(pz,2)==1);
+if ergodic && all(isfinite(pz))
 
-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;
+	TlogT=T.*slog(T);
+	HXZ=-sum(TlogT,1)*pz;
 
-IXZ=max(0,(sum(TlogT)-slog(pz)')*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;
 
-Info=[HXZ;IXZ;IXYZ];
-if any(~isreal(Info))
-	disp('unreal information');
-	keyboard
-elseif any(Info<0)
-	if any(Info<-eps)
-		disp('negative information');
-		disp(Info);
-		keyboard
-	else
-		Info=max(0,Info);
+	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)

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/private/mc_global_info1.m	Sun Feb 26 23:11:10 2012 +0000
@@ -0,0 +1,82 @@
+function [Info,ergodic]=mc_global_info1(T,tol)
+% 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.
+
+
+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
+
+function y=slog(x)
+% slog - safe log, x<=0 => slog(x)=0
+%
+% slog :: real -> real.
+
+	x(x<=0)=1;
+	y=log(x);
+
+% 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

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/private/randnat.m	Sun Feb 26 23:11:10 2012 +0000
@@ -0,0 +1,18 @@
+function X=randnat(M,varargin)
+% randnat - Draw Random natural numbers from uniform distribution
+%
+% randnat :: 
+%    M:natural   ~'maximum value', 
+%    D:[[1,E]]   ~'size of E-dimensional array'
+% -> [[D]->[M]] ~'size D array of numbers in 1..M'.
+%
+% The size is specified as in the RAND function, except that if
+% only one dimension is given, we create a vector, not a square
+% matrix, eg
+%
+% randnat(10,3) :: [[3,3]->[10]].
+% randnat(20,2,5,7) :: [[2,5,7]->[20]].
+% randnat(20,5) :: [[5,1]->[20]].
+
+X=1+floor(M*rand(tosize(varargin{:})));
+

--- a/private/sample_transmat_hdp.m	Wed Feb 22 21:24:44 2012 +0000
+++ b/private/sample_transmat_hdp.m	Sun Feb 26 23:11:10 2012 +0000
@@ -1,4 +1,4 @@
-function [T,I]=sample_transmat_hdp(A0,B0,A1,B1,K,L,T0)
+function T=sample_transmat_hdp(A0,B0,A1,B1,K,L,T0)
 % sample_transmat_hdp - Sample 1st order Markov transition matrix using HDP
 %
 % sample_transmat_hdp ::
@@ -8,8 +8,7 @@
 %    nonneg ~'Scale parameter for sub Dirichlet Gamma prior',
 %    natural ~'size of alphabet',
 %    natural ~'num of transition matrices to produce',
-% -> [[K,K,L]] ~ 'lots of transition matrices',
-%    [[L,3]]   ~ 'entropy rate, redundancy, and predinfo rate for all'.
+% -> [[K,K,L]] ~ 'lots of transition matrices'.
 %
 % An extra optional argument of type [[K,K]] can be added, which
 % is used as a sort of bias for the sampled matrices.
@@ -26,14 +25,7 @@
 	T=sample_dir(repmat(Alpha1,[K,K,1]).*HT);
 
 	if any(sum(H0)<0.99)
-		disp('normalisation problem in H')
-		keyboard
-	end
-	if nargout>=2
-		I=zeros(K,3);
-		for i=1:L
-			I(i,:)=mc_global_info(T(:,:,i));
-		end
+		error('normalisation problem in H')
 	end
 end