def randomwrap(func):
    """
    Decorator for random value generators

    Allows passing of sequence of parameters, as well as a size argument.


    Convention:

      - If size=1 and the parameters are all scalars, return a scalar.
      - If size=1, the random variates are 1D.
      - If the parameters are scalars and size > 1, the random variates are 1D.
      - If size > 1 and the parameters are sequences, the random variates are
        aligned as (size, max(length)), where length is the parameters size.


    :Example:
      >>> rbernoulli(.1)
      0
      >>> rbernoulli([.1,.9])
      np.asarray([0, 1])
      >>> rbernoulli(.9, size=2)
      np.asarray([1, 1])
      >>> rbernoulli([.1,.9], 2)
      np.asarray([[0, 1],
             [0, 1]])
    """

    # Find the order of the arguments.
    refargs, defaults = utils.get_signature(func)
    # vfunc = np.vectorize(self.func)
    npos = len(refargs) - len(defaults)  # Number of pos. arg.
    nkwds = len(defaults)  # Number of kwds args.
    mv = func.__name__[
        1:] in mv_continuous_distributions + mv_discrete_distributions

    # Use the NumPy random function directly if this is not a multivariate
    # distribution
    if not mv:
        return func

    def wrapper(*args, **kwds):
        # First transform keyword arguments into positional arguments.
        n = len(args)
        if nkwds > 0:
            args = list(args)
            for i, k in enumerate(refargs[n:]):
                if k in kwds.keys():
                    args.append(kwds[k])
                else:
                    args.append(defaults[n - npos + i])

        r = []
        s = []
        largs = []
        nr = args[-1]
        length = [np.atleast_1d(a).shape[0] for a in args]
        dimension = [np.atleast_1d(a).ndim for a in args]
        N = max(length)
        if len(set(dimension)) > 2:
            raise('Dimensions do not agree.')
        # Make sure all elements are iterable and have consistent lengths, ie
        # 1 or n, but not m and n.

        for arg, s in zip(args, length):
            t = type(arg)
            arr = np.empty(N, type)
            if s == 1:
                arr.fill(arg)
            elif s == N:
                arr = np.asarray(arg)
            else:
                raise RuntimeError('Arguments size not allowed: %s.' % s)
            largs.append(arr)

        if mv and N > 1 and max(dimension) > 1 and nr > 1:
            raise ValueError(
                'Multivariate distributions cannot take s>1 and multiple values.')

        if mv:
            for i, arg in enumerate(largs[:-1]):
                largs[0] = np.atleast_2d(arg)

        for arg in zip(*largs):
            r.append(func(*arg))

        size = arg[-1]
        vec_stochastics = len(r) > 1
        if mv:
            if nr == 1:
                return r[0]
            else:
                return np.vstack(r)
        else:
            if size > 1 and vec_stochastics:
                return np.atleast_2d(r).T
            elif vec_stochastics or size > 1:
                return np.concatenate(r)
            else:  # Scalar case
                return r[0][0]

    wrapper.__doc__ = func.__doc__
    wrapper.__name__ = func.__name__
    return wrapper
Decorator for random value generators

    Allows passing of sequence of parameters, as well as a size argument.


    Convention:

      - If size=1 and the parameters are all scalars, return a scalar.
      - If size=1, the random variates are 1D.
      - If the parameters are scalars and size > 1, the random variates are 1D.
      - If size > 1 and the parameters are sequences, the random variates are
        aligned as (size, max(length)), where length is the parameters size.


    :Example:
      >>> rbernoulli(.1)
      0
      >>> rbernoulli([.1,.9])
      np.asarray([0, 1])
      >>> rbernoulli(.9, size=2)
      np.asarray([1, 1])
      >>> rbernoulli([.1,.9], 2)
      np.asarray([[0, 1],
             [0, 1]])