code: def rdirichlet(theta, size=1):
      
      gammas = np.vstack([rgamma(theta, 1) for i in xrange(size)])
      if size > 1 and np.size(theta) > 1:
          return (gammas.T / gammas.sum(1))[:-1].T
      elif np.size(theta) > 1:
          return (gammas[0] / gammas[0].sum())[:-1]
      else:
          return 1.
docs: Dirichlet random variates.