comparison private/sample_transmat_hdp.m @ 0:be936975f254

Initial check in.

+function [T,I]=sample_transmat_hdp(A0,B0,A1,B1,K,L,T0)
+% sample_transmat_hdp - Sample 1st order Markov transition matrix using HDP
+%
+% sample_transmat_hdp ::
+%    nonneg ~'Shape parameter for base Dirichlet Gamma prior',
+%    nonneg ~'Scale parameter for base Dirichlet Gamma prior',
+%    nonneg ~'Shape parameter for sub Dirichlet Gamma prior',
+%    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'.
+%
+% An extra optional argument of type [[K,K]] can be added, which
+% is used as a sort of bias for the sampled matrices.
+
+	if nargin<7, T0=ones(K,K); end
+
+	Alpha0=gamrnd(A0,B0,1,1,L);
+	Alpha1=gamrnd(A1,B1,1,1,L);
+
+
+	% base distributions
+	H0=sample_dir(repmat(Alpha0/K,[K,1,1]));
+	HT=stoch(repmat(H0,[1,K,1])).*repmat(T0/sum(T0(:)),[1,1,L]);
+	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
+	end
+end
+
+% taken from stats/models/dirichlet.m>stoch_gamma
+function y=sample_dir(beta)
+	x=randg(beta);
+	z=sum(x,1);
+	y=x./repmat(z,size(x,1),1);
+	bad=find(z==0);
+	for i=1:length(bad)
+		beta_i=beta(:,bad(i)); % misbehaving params
+		k=(beta_i==max(beta_i)); % entries with max beta
+		y(:,bad(i))=mnrnd(1,stoch(k));
+	end
+end
+
+function [Y,Z]=stoch(X),
+	Z=sum(X,1);
+	Y=X./repmat(Z,size(X,1),1);
+end
+
+