#GENERATE SAMPLE FROM DIRICHLET--------------------------
rdirichlet _ function ( n, p ) {
# return n random samples from a Dirichlet distribution with parameter p
  if (   !is.vector(n, "numeric")
       | length(n) != 1
       | !is.vector(p, "numeric")
     ) { stop("error in call to rdirichlet") }
  mat _ matrix ( NA, n, length(p) )
  mat[,1] _ rbeta ( n, p[1], sum(p[-1]) )
  
for ( i in 2:(length(p)-1) ) { 
 
  mat[,i] _  (    rbeta ( n, p[i], sum(p[(i+1):length(p)]) )
                  * ( 1 - apply ((mat[,1:(i-1),drop=F]), 1, sum) )
               )
  }
  mat[,length(p)] _  1 - apply ( (mat[,-length(p),drop=F]), 1, sum )

  return ( mat )
}