This is mostly about recording the conversion between Numpy/Scipy
    conventions and Wikipedia conventions. Some equations:

    mean = alpha / beta
    variance = alpha / beta**2
    mode = (alpha - 1) / beta [if alpha > 1; otherwise undefined]
    skewness = 2 / sqrt(alpha)

def sample_gamma(alpha, beta, size):
    """This is mostly about recording the conversion between Numpy/Scipy
    conventions and Wikipedia conventions. Some equations:

    mean = alpha / beta
    variance = alpha / beta**2
    mode = (alpha - 1) / beta [if alpha > 1; otherwise undefined]
    skewness = 2 / sqrt(alpha)
    """

    if alpha <= 0:
        raise ValueError('alpha must be positive; got %e' % alpha)
    if beta <= 0:
        raise ValueError('beta must be positive; got %e' % beta)
    return np.random.gamma(alpha, scale=1./beta, size=size)