rdirichlet <- function( n, alpha ) {

  p <- length( alpha )

  g <- matrix( NA, n, p )

  for ( i in 1:n ) {

    for ( j in 1:p ) {

      g[ i, j ] <- rgamma( 1, rate = 1, shape = alpha[ j ] )

    }

    g[ i, ] <- g[ i, ] / sum( g[ i, ] )

  }

  return( g )

}