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