wolffd@0: function x = dirichletrnd(alpha)
wolffd@0: %DIRICHLETRND Random vector from a dirichlet distribution.
wolffd@0: %   x = dirichletrnd(alpha) returns a vector randomly selected
wolffd@0: %   from the Dirichlet distribution with parameter vector alpha.
wolffd@0: %
wolffd@0: %   The algorithm used is the following:
wolffd@0: %   For each alpha(i), generate a value s(i) with distribution
wolffd@0: %   Gamma(alpha(i),1).  Now x(i) = s(i) / sum_j s(j).
wolffd@0: %   
wolffd@0: %   The above algorithm was recounted to me by Radford Neal, but
wolffd@0: %   a reference would be appreciated...
wolffd@0: %   Do the gamma parameters always have to be 1?
wolffd@0: %
wolffd@0: %   Author: David Ross
wolffd@0: %   $Id: dirichletrnd.m,v 1.1.1.1 2005/05/22 23:32:12 yozhik Exp $
wolffd@0: 
wolffd@0: %-------------------------------------------------
wolffd@0: % Check the input
wolffd@0: %-------------------------------------------------
wolffd@0: error(nargchk(1,1,nargin));
wolffd@0: 
wolffd@0: if min(size(alpha)) ~= 1 | length(alpha) < 2
wolffd@0:     error('alpha must be a vector of length at least 2');
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: %-------------------------------------------------
wolffd@0: % Main
wolffd@0: %-------------------------------------------------
wolffd@0: gamma_vals = gamrnd(alpha, ones(size(alpha)), size(alpha));
wolffd@0: denom = sum(gamma_vals);
wolffd@0: x = gamma_vals / denom;