Random beta-binomial variates.

def rbetabin(alpha, beta, n, size=None):
    """
    Random beta-binomial variates.
    """

    phi = np.random.beta(alpha, beta, size)
    return np.random.binomial(n, phi)