def DiscreteGamma(alpha, beta, ncats):
alpha = float(alpha)<EOL>beta = float(beta)<EOL>assert alpha > <NUM_LIT:0><EOL>assert beta > <NUM_LIT:0><EOL>assert ncats > <NUM_LIT:1><EOL>scale = <NUM_LIT:1.0> / beta<EOL>catmeans = scipy.ndarray(ncats, dtype='<STR_LIT:float>')<EOL>for k in range(ncats):<EOL><INDENT>if k == <NUM_LIT:0>:<EOL><INDENT>lower = <NUM_LIT:0.0><EOL>gammainc_lower = <NUM_LIT:0.0><EOL><DEDENT>else:<EOL><INDENT>lower = upper<EOL>gammainc_lower = gammainc_upper<EOL><DEDENT>if k == ncats - <NUM_LIT:1>:<EOL><INDENT>upper = float('<STR_LIT>')<EOL>gammainc_upper = <NUM_LIT:1.0><EOL><DEDENT>else:<EOL><INDENT>upper = scipy.stats.gamma.ppf((k + <NUM_LIT:1>) / float(ncats), alpha,<EOL>scale=scale)<EOL>gammainc_upper = scipy.special.gammainc(alpha + <NUM_LIT:1>, upper * beta)<EOL><DEDENT>catmeans[k] = alpha * ncats * (gammainc_upper - gammainc_lower) / beta<EOL><DEDENT>assert scipy.allclose(catmeans.sum() / ncats, alpha / beta), (<EOL>"<STR_LIT>"<EOL>"<STR_LIT>").format(catmeans,<EOL>catmeans.sum() / ncats, alpha, beta, alpha / beta)<EOL>return catmeans<EOL>
Returns category means for discretized gamma distribution.

    The distribution is evenly divided into categories, and the
    mean of each category is returned.

    Args:
        `alpha` (`float` > 0)
            Shape parameter of gamma distribution.
        `beta` (`float` > 0)
            Inverse scale parameter of gamma distribution.
        `ncats` (`int` > 1)
            Number of categories in discretized gamma distribution.

    Returns:
        `catmeans` (`scipy.ndarray` of `float`, shape `(ncats,)`)
            `catmeans[k]` is the mean of category `k` where
            `0 <= k < ncats`.

    Check that we get values in Figure 1 of Yang, J Mol Evol, 39:306-314
    >>> catmeans = DiscreteGamma(0.5, 0.5, 4)
    >>> scipy.allclose(catmeans, scipy.array([0.0334, 0.2519, 0.8203, 2.8944]), atol=1e-4)
    True

    Make sure we get expected mean of alpha / beta
    >>> alpha = 0.6
    >>> beta = 0.8
    >>> ncats = 6
    >>> catmeans = DiscreteGamma(alpha, beta, ncats)
    >>> scipy.allclose(catmeans.sum() / ncats, alpha / beta)
    True