function r = drchrnd(a,n)
%{
take a sample from a dirichlet distribution
copied from Mathworks website with a few changes

note that length(a) vectors are returned, you probably don't want the last
because they add to one.
Dick Startz
March 2017

%}
p = length(a);
if nargin == 1
    r = gamrnd(a,1);
    r = r ./sum(r,2);
else
    r = gamrnd(repmat(a,n,1),1,n,p);
    r = r ./ repmat(sum(r,2),1,p);
end
end
