view 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