rdirichlet <- function( m, alpha ) {

  k <- length( alpha )

  theta.star <- matrix( 0, m, k )

  for ( j in 1:k ) {

    theta.star[ , j ] <- rgamma( m, alpha[ j ], 1 )

  }

  theta.star <- theta.star / apply( theta.star, 1, sum )

  return( theta.star )

}