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