Численная стабильность AB-тестирования в Python

Question

Учитывая X~Beta(a1, b1) и Y~Beta(a2, b2) Я хочу рассчитать E[f(x, y)], для AB-тестирования в python численно стабильным способом.

Чтобы получить первое впечатление, я использовал f(x, y) = I(x>y) и ожидаемые значения(приблизительно) 0,5 для симметричных распределений, то есть a1=b1 и a2=b2.Код, кажется, работает для малых значений, но разбивается на более крупные.

Вопросы, таким образом:

1. Почему приведенный ниже код работает только для небольших значений?
2. Как будет выглядеть численно более стабильное решение?

Я также пробовал версии плотности, не работающие с журналами, и прямые решения, такие как это , но никогда не достигал численной стабильности.

from scipy.special import betaln
from math import exp, log
from scipy.integrate import dblquad

def loss(x, y):
    return 1 if x>y else 0


def density(x, y, a1, b1, a2, b2):
    return exp((a1-1) * log(x) + (b1-1) * log(1-x) + (a2-1) * log(y) + (b2-1) * log(1-y) - lbeta(a1, b1) - lbeta(a2, b2))


def expected_loss_dblquad(a1, b1, a2, b2, *args, **kwargs):
    return dblquad(lambda x, y: loss(x, y) * density(x, y, a1, b1, a2, b2),
                   0, 1, lambda x: 0, lambda x: 1, *args, **kwargs)

expected_loss_dblquad(1, 1, 1, 1)  # returns (0.499997642156207, 1.4849469934935522e-08)

expected_loss_dblquad(10, 10, 100, 100)  # returns (0.17088036863463973, 8.33972164867105e-08)