Эффективная нижняя неполная бета-функция

Кто-нибудь реализовал эффективную по времени и относительно точную функцию нижней неполной бета-версии? Моя текущая попытка дает несколько неточные результаты, и если я уменьшу значение приращения, время ожидания сценария истечет; включение функции «как есть» в другие сценарии приводит к их тайм-ауту из-за того, что для каждой точки данных рассчитывается длинный цикл.

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

Суть в том, чтобы использовать эту функцию для вычисления точных значений P в pinescript с использованием распределения Стьюдента, которое требует CDF бета, который, в свою очередь, требует этой функции.

_beta_incomplete(x, a, b) =>
    var float i = 0.000001
    float _a = a - 1
    float _b = b - 1
    float integral = 0.0
    for t = i to x by i
        integral := integral + pow(t, _a) * pow(1 - t, _b)
    integral * i

ОБНОВЛЕНИЕ: Моя попытка воспроизвести эту реализацию. Вывод просто неправильный, но я не уверен, где я ошибся.

// The functions gamma_incomplete, standard_normal_pdf and standard_normal_cdf have already been vetted
_h(z, x2, c, v) =>
    float q = v * pow(z + c, 2) / x2
    -log(2) - 0.5 * (q - v - (v - 2) * log(q / v) + log(v) + log(2 * pi) + pow(z, 2))

_g(z, x2, c, v) => _gamma_incomplete(v * pow(z + c, 2) / (2 * x2), v / 2) * standard_normal_pdf(z)

_t_distribution_cdf(x, dof) =>
    var float er = pow(10, -16)
    var float ler = log(er)
    var float th = 37.5194
    // Non-centrality parameter...don't know what it's for, so set to 0
    var float c = 0
    float x2 = pow(x, 2)
    float sdof = sqrt(dof)
    float zmod = (-c * (x2 + 2 * dof) + x * sqrt(4 * dof * (dof - 2) + x2 * (pow(c, 2) + 4 * (dof - 2)))) / (2 * (x2 + dof))
    float qer = dof + 2 * er + 1.62 * sqrt(dof * er) + 0.63012 * sdof * ler - 1.12032 * sdof - 2.48 * sqrt(er) - 0.65381 * ler - 0.22872
    // TODO: Is "ea" the right term to use here?
    float ea = exp(_h(zmod, x2, c, dof) + ler)
    float a = max(max(max(-c, -th), sqrt(x2 * abs(qer) / dof) - c), ea)
    float b = min(th, ea)
    // This is my attempt to implement "nsubs" (page 7). I think they want to
    // "bin" the integrands into 16 separate calculations which are summed to the integral
    float factor = abs(a - b) / 16
    float integral = standard_normal_cdf(a)
    for z = a to b by factor
        integral := integral + _g(z, x2, c, dof)
    integral