def multinomial_helper(random_state, n, pvals, size):
    # 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 = numpy.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 = numpy.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 = numpy.sum(pvi)
        elif pvi[-1] < 5e-5:  # will this even work?
            pvi = pvi * (1.0 - 5e-5)
            pisum = numpy.sum(pvi)
        assert pisum <= 1.0, pisum
        out[mi] = random_state.multinomial(n=n[ni],
                                           pvals=pvi.astype('float64'))
    return out
Helper function drawing from multinomial distributions.
dataset/ETHPy150Open rizar/attention-lvcsr/libs/Theano/theano/tensor/raw_random.py/multinomial_helper