code: def multinomial_helper(random_state, n, pvals, size):
      """
      Helper function drawing from multinomial distributions.
  
      This is a generalization of numpy.random.multinomial to the case where
      n and pvals are tensors.
  
      """
      # Figure out the shape if it's None
      # Note: the output ndim will be ndim+1, because the multinomial
      # adds a dimension. The length of that dimension is pvals.shape[-1].
      if size is not None:
          ndim = len(size)
      else:
          ndim = max(n.ndim, pvals.ndim - 1)
  
      # broadcast n to ndim dimensions and pvals to ndim+1
      if n.ndim > ndim:
          raise ValueError(
              "n.ndim (%i) should not be larger than len(size) (%i)"
              % (n.ndim, ndim),
              n,
              size,
          )
      if n.ndim < ndim:
          n = n.reshape((1,) * (ndim - n.ndim) + n.shape)
  
      if pvals.ndim - 1 > ndim:
          raise ValueError(
              "pvals.ndim-1 (%i) should not be larger than len(size) (%i)"
              % (pvals.ndim - 1, ndim),
              pvals,
              size,
          )
      if pvals.ndim - 1 < ndim:
          pvals = pvals.reshape(
              (1,) * (ndim - pvals.ndim + 1) + pvals.shape
          )
  
      if size is not None:
          size = tuple(size)
      else:
          size = ()
          for dim in xrange(ndim):
              dim_len = max(n.shape[dim], pvals.shape[dim])
              size = size + (dim_len,)
      out_size = size + (pvals.shape[-1],)
  
      # Build the indices over which to loop
      # Note that here, the rows (inner-most 1D subtensors) of pvals and out
      # are indexed, not their individual elements
      out = np.ndarray(out_size)
      broadcast_ind = _generate_broadcasting_indices(
          size, n.shape, pvals.shape[:-1]
      )
      # Iterate over these indices, drawing from one multinomial at a
      # time from numpy
      assert pvals.min() >= 0
      for mi, ni, pi in zip(*broadcast_ind):
          pvi = pvals[pi]
  
          # This might someday be fixed upstream
          # Currently numpy raises an exception in this method if the sum
          # of probabilities meets or exceeds 1.0.
          # In  perfect arithmetic this would be correct, but in float32 or
          # float64 it is too strict.
          pisum = np.sum(pvi)
          if (
              1.0 < pisum < 1.0 + 1e-5
          ):  # correct if we went a little over
              # because mtrand.pyx has a ValueError that will trigger if
              # sum(pvals[:-1]) > 1.0
              pvi = pvi * (1.0 - 5e-5)
              # pvi = pvi * .9
              pisum = np.sum(pvi)
          elif pvi[-1] < 5e-5:  # will this even work?
              pvi = pvi * (1.0 - 5e-5)
              pisum = np.sum(pvi)
          assert pisum <= 1.0, pisum
          out[mi] = random_state.multinomial(
              n=n[ni], pvals=pvi.astype("float64")
          )
      return out
  
imports: from __future__ import absolute_import, print_function, division
  import sys
  from copy import copy
  import numpy
  from six import string_types
  from six.moves import reduce, xrange
  import theano
  from theano import tensor
  from theano.tensor import opt
  from theano import gof
  from theano.compile import optdb
label_desc: Use image node_id 13 for calling a global function with example usage: multinomial_helper(random_state, n, pvals, size) and returns: out