
function r = dirichlet_sample_good(a,n)
  if nargin < 2
    n = 1;
  end
  row = (size(a,1) == 1);
  a = a(:);
  y = gamrnd(repmat(a, 1, n),1);
  % randgamma is faster
  %y = randgamma(repmat(a, 1, n));
  r = sum(y,1);
  r(find(r == 0)) = 1;
  r = y./repmat(r, size(y,1), 1);
  if row
    r = r';
  end  
