%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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

