Has anyone implemented a time-efficient and relatively accurate lower incomplete beta function? My current attempt yields somewhat inaccurate results and if I lower the increment value the script times out; including the function as-is in other scripts times them out due to the long loop being calculated on every data point.
Looking at solutions to other somewhat related problems I can't help but feel there's a similar tact to take with this but I don't understand why their solutions work (using seemingly arbitrary floating point values). An ideal solution would be to minimize the use of loops somehow.
The longshot here is to use this function for the sake of calculating accurate P values in pinescript using the student's t distribution which requires the CDF of beta which in turns requires this function.
_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
UPDATE: My attempt to replicate this implementation. The output is simply wrong, but I'm not sure where I've gone wrong.
// 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
