wolffd@0: function theta = dirichlet_sample(alpha, N)
wolffd@0: % SAMPLE_DIRICHLET Sample N vectors from Dir(alpha(1), ..., alpha(k))
wolffd@0: % theta = sample_dirichlet(alpha, N)
wolffd@0: % theta(i,j) = i'th sample of theta_j, where theta ~ Dir
wolffd@0: 
wolffd@0: % We use the method from p. 482 of "Bayesian Data Analysis", Gelman et al.
wolffd@0: 
wolffd@0: assert(alpha > 0);
wolffd@0: k = length(alpha);
wolffd@0: theta = zeros(N, k);
wolffd@0: scale = 1; % arbitrary
wolffd@0: for i=1:k
wolffd@0:   %theta(:,i) = gamrnd(alpha(i), scale, N, 1);
wolffd@0:   theta(:,i) = gamma_sample(alpha(i), scale, N, 1);
wolffd@0: end
wolffd@0: %theta = mk_stochastic(theta);
wolffd@0: S = sum(theta,2); 
wolffd@0: theta = theta ./ repmat(S, 1, k);