To Generate Random Numbers from a Dirichlet Distribution

The following code snippet is copied from the MATLAB Topic Modeling Toolbox by Mark Steyvers and Tom Griffiths:

function r = drchrnd(a,n)
% take a sample from a dirichlet distribution
p = length(a);
r = gamrnd(repmat(a,n,1),1,n,p);
r = r ./ repmat(sum(r,2),1,p);

The following is an example that generates three discrete distributions from a symmetric Dirichlet distribution Dir( \theta ; [ 1 1 1 1 ] ):

>> A = drchrnd([1 1 1 1], 3)

A =

0.3889 0.1738 0.0866 0.3507
0.0130 0.0874 0.6416 0.2579
0.0251 0.0105 0.2716 0.6928

>> sum(A, 2)

ans =

1
1
1