#----------------------Subroutines-------------------------------------------------
#The function rmultinom returns a random multinomial vector
#n=number of variables desired
#size=sum of cells
#prob=vector of cell probabilities
rmultinom <- function(n, size, prob) {
    K <- length(prob) # #{classes}
    matrix(tabulate(sample(K, n * size, repl = TRUE, prob) + K *0:(n - 1),
                    nbins = n * K),
           nrow = n, ncol = K, byrow = TRUE)}
#----------------------------------------------------------------------------------
rdirichlet <- function(n, shape) {
  d <- length(shape)
  x <- rgamma(n*d, rep(shape, n))
  x <- matrix(x, nrow = d)
  sweep(x, 2, apply(x, 2, sum), "/")
}

#So for example, to get 1000 Dirichlets with common parameter vector
#c(5,10,2,3), you type
#
#rdirichlet(1000, c(5,10,2,3))
#
#This returns a 4 by 1000 matrix, with each column representing a
#single Dirichlet observation.
#----------------------------------------------------------------------------------