%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Draws form the Dirichlet
%%%%%%%%%%%%%%%%%%%%%%%%%%
function D=dirrnd(a,R);
%
% Draws of random vector(s) from the Dirichlet distribution
% with parameter vector "a (k by 1)"
% R denotes the number of desired draws
% D is a k by R matrix, with k=length(a)
% I'm using the more efficient approach suggested in Gelman et al, p.
% 582, based on Gamma draws.
% This method performed better than the alternative based on Beta draws -
% see script 'dirichlet_draws_test'
k=length(a);
D=zeros(k,R); %R draws of Dirichlet vectors
for i=1:R
int=gamrnd(a,1); %returns k by 1 vector of gamma draws with
%pairwise-common shape & scale
D(:,i)=int/sum(int);
end