code: def rmultinomial(n, p, size=None):
      
      # Leaving size=None as the default means return value is 1d array
      # if not specified-- nicer.
  
      # Single value for p:
      if len(np.shape(p)) == 1:
          return np.random.multinomial(n, p, size)
  
      # Multiple values for p:
      if np.isscalar(n):
          n = n * np.ones(np.shape(p)[0], dtype=np.int)
      out = np.empty(np.shape(p))
      for i in xrange(np.shape(p)[0]):
          out[i, :] = np.random.multinomial(n[i], p[i,:], size)
      return out
docs: Random multinomial variates.