Quite fast and accurate:

    from scipy.stats import beta
    import numpy as np

    def betamax(params0, params1, integration_points=5000):
        u0 = params0[0] / sum(params0)
        u1 = params1[0] / sum(params1)
        flip = u0 < u1
        if flip:
            params0, params1 = params1, params0
        q = beta(*params0).ppf(np.linspace(0, 1, integration_points + 2)[1:-1])
        c = beta(*params1).cdf(q)
        p = c.mean()
        if flip:
            return 1 - p
        return p