Repository: scipy/scipy
File: scipy/stats/_continuous_distns.py (Python)

Code before:
            t3 = -50/b - 10*b**-2.0 - b**-3.0 + b**-4.0
            return log_term + (t1 + t2 + t3) / 120

        def asymptotic_b_large(a, b):
            sum_ab = a + b
            t1 = sc.gammaln(a) - (a - 1) * sc.psi(a)
            t2 = (
                - 1/(2*b) + 1/(12*b) - b**-2.0/12 - b**-3.0/120 + b**-4.0/120
                + b**-5.0/252 - b**-6.0/252 + 1/sum_ab - 1/(12*sum_ab)
                + sum_ab**-2.0/6 + sum_ab**-3.0/120 - sum_ab**-4.0/60
                - sum_ab**-5.0/252 + sum_ab**-6.0/126
            )
            log_term = sum_ab*np.log1p(a/b) + np.log(b) - 2*np.log(sum_ab)
            return t1 + t2 + log_term

        def asymptotic_a_large(a, b):
            return asymptotic_b_large(b, a)

        def threshold_large(v):
            j = np.floor(np.log10(v))
            d = np.floor(v / 10 ** j) + 2
            return xpx.apply_where(v != 1.0, (d, j), lambda d_, j_: d_ * 10**(7 + j_),
                                   fill_value=1000)

        threshold_a = threshold_large(a)
        threshold_b = threshold_large(b)
        return _lazyselect([np.logical_and(a >= 4.96e6, b >= 4.96e6),
                            np.logical_and(np.logical_and(a <= 4.9e6, b - a >= 1e6),
                                           b >= threshold_a),
                            np.logical_and(np.logical_and(b <= 4.9e6, a - b >= 1e6),
                                           a >= threshold_b),
                            np.logical_and(a < 4.9e6, b < 4.9e6)
                           ],
                           [asymptotic_ab_large, asymptotic_b_large,
                            asymptotic_a_large, regular],
                           [a, b]
        )


beta = beta_gen(a=0.0, b=1.0, name='beta')


class betaprime_gen(rv_continuous):
    r"""A beta prime continuous random variable.

    %(before_notes)s

    Notes
    -----
    The probability density function for `betaprime` is:


Review comment (suggestion):
I don't think it matters much, but it looks like the most direct translation would have been to vectorize `threshold_large` and just use `apply_where` (twice) to assign `threshold_a` and `threshold_b`, right?
This is fine, it just might be nice to use `threshold_large` as the callable of the `apply_where`.

Code after:
            t3 = -50/b - 10*b**-2.0 - b**-3.0 + b**-4.0
            return log_term + (t1 + t2 + t3) / 120

        def asymptotic_b_large(a, b):
            sum_ab = a + b
            t1 = sc.gammaln(a) - (a - 1) * sc.psi(a)
            t2 = (
                - 1/(2*b) + 1/(12*b) - b**-2.0/12 - b**-3.0/120 + b**-4.0/120
                + b**-5.0/252 - b**-6.0/252 + 1/sum_ab - 1/(12*sum_ab)
                + sum_ab**-2.0/6 + sum_ab**-3.0/120 - sum_ab**-4.0/60
                - sum_ab**-5.0/252 + sum_ab**-6.0/126
            )
            log_term = sum_ab*np.log1p(a/b) + np.log(b) - 2*np.log(sum_ab)
            return t1 + t2 + log_term

        def asymptotic_a_large(a, b):
            return asymptotic_b_large(b, a)

        def threshold_large(v):
            j = np.floor(np.log10(v))
            d = np.floor(v / 10 ** j) + 2
            return xpx.apply_where(v != 1.0, (d, j), lambda d_, j_: d_ * 10**(7 + j_),
                                   fill_value=1000)

        threshold_a = threshold_large(a)
        threshold_b = threshold_large(b)
        return _lazyselect([(a >= 4.96e6) & (b >= 4.96e6),
                            (a <= 4.9e6) & (b - a >= 1e6) & (b >= threshold_a),
                            (b <= 4.9e6) & (a - b >= 1e6) & (a >= threshold_b),
                            (a < 4.9e6) & (b < 4.9e6)
                           ],
                           [asymptotic_ab_large, asymptotic_b_large,
                            asymptotic_a_large, regular],
                           [a, b]
        )


beta = beta_gen(a=0.0, b=1.0, name='beta')


class betaprime_gen(rv_continuous):
    r"""A beta prime continuous random variable.

    %(before_notes)s

    Notes
    -----
    The probability density function for `betaprime` is:

    .. math::