Files

128 lines
4.6 KiB
Plaintext

// Licensed under the Apache License, Version 2.0
// © mihakralj
//@version=6
indicator("F-Distribution CDF (FDIST)", "FDIST", overlay=false, precision=6)
//@function Log-gamma via Lanczos approximation (g=7, 9 coefficients)
//@param z Input value (z > 0)
//@returns ln(Gamma(z))
lnGamma(simple float z) =>
float g = 7.0
array<float> c = array.from(
0.99999999999980993,
676.5203681218851,
-1259.1392167224028,
771.32342877765313,
-176.61502916214059,
12.507343278686905,
-0.13857109526572012,
9.9843695780195716e-6,
1.5056327351493116e-7)
float zz = z - 1.0
float x = array.get(c, 0)
for i = 1 to 8
x += array.get(c, i) / (zz + i)
float t = zz + g + 0.5
0.5 * math.log(2.0 * math.pi) + (zz + 0.5) * math.log(t) - t + math.log(x)
//@function Regularized incomplete beta function I_x(a,b) via continued fraction (Lentz)
//@param x Evaluation point (0 <= x <= 1)
//@param a Shape parameter alpha (a > 0)
//@param b Shape parameter beta (b > 0)
//@returns CDF value P(X <= x) for Beta(a,b)
betaReg(series float x, simple float a, simple float b) =>
if x <= 0.0
0.0
else if x >= 1.0
1.0
else
float lnPfx = a * math.log(x) + b * math.log(1.0 - x) - math.log(a) - lnGamma(a) - lnGamma(b) + lnGamma(a + b)
float front = math.exp(lnPfx)
bool flip = x > (a + 1.0) / (a + b + 2.0)
float xx = flip ? 1.0 - x : x
float aa = flip ? b : a
float bb = flip ? a : b
float lnPfx2 = aa * math.log(xx) + bb * math.log(1.0 - xx) - math.log(aa) - lnGamma(aa) - lnGamma(bb) + lnGamma(aa + bb)
float front2 = math.exp(lnPfx2)
float TINY = 1e-30
float EPS = 1e-10
int MAXITER = 200
float f = TINY
float C = TINY
float D = 0.0
float delta = 0.0
for m = 0 to MAXITER - 1
float d_val = 0.0
int mm = m / 2
if m == 0
d_val := 1.0
else if m % 2 == 0
float mf = float(mm)
d_val := mf * (bb - mf) * xx / ((aa + 2.0 * mf - 1.0) * (aa + 2.0 * mf))
else
float mf = float(mm) + 1.0
d_val := -(aa + mf - 1.0) * (aa + bb + mf - 1.0) * xx / ((aa + 2.0 * mf - 2.0) * (aa + 2.0 * mf - 1.0))
D := 1.0 + d_val * D
if math.abs(D) < TINY
D := TINY
D := 1.0 / D
C := 1.0 + d_val / C
if math.abs(C) < TINY
C := TINY
delta := C * D
f *= delta
if math.abs(delta - 1.0) < EPS
break
float result = front2 * f
flip ? 1.0 - result : result
//@function Computes F-Distribution CDF for a normalized price series
//@param source Series to transform
//@param period Lookback period for min-max normalization
//@param d1 Numerator degrees of freedom (d1 > 0)
//@param d2 Denominator degrees of freedom (d2 > 0)
//@returns CDF value in [0,1]: P(F <= x) via regularized incomplete beta
//@optimized Lentz continued fraction converges in ~10-20 iterations
fdist(series float source, simple int period, simple float d1, simple float d2) =>
if period <= 0
runtime.error("Period must be greater than 0")
if d1 <= 0.0
runtime.error("d1 must be greater than 0")
if d2 <= 0.0
runtime.error("d2 must be greater than 0")
float minVal = source
float maxVal = source
for i = 1 to period - 1
float v = source[i]
if not na(v)
if v < minVal
minVal := v
if v > maxVal
maxVal := v
float range = maxVal - minVal
float x = range > 0.0 ? (source - minVal) / range : 0.5
float safeX = math.max(0.0, x)
float t = d1 * safeX / (d1 * safeX + d2)
betaReg(t, d1 / 2.0, d2 / 2.0)
// ---------- Main loop ----------
// Inputs
i_source = input.source(close, "Source")
i_period = input.int(50, "Lookback Period", minval=2, maxval=5000, tooltip="Min-max normalization window")
i_d1 = input.float(5.0, "d1 (numerator df)", minval=0.1, step=1.0, tooltip="Numerator degrees of freedom")
i_d2 = input.float(5.0, "d2 (denominator df)", minval=0.1, step=1.0, tooltip="Denominator degrees of freedom")
// Calculation
float result = fdist(i_source, i_period, i_d1, i_d2)
// Plot
plot(result, "FDIST", color=color.yellow, linewidth=2)
hline(0.5, "Midline", color=color.gray, linestyle=hline.style_dotted)
hline(0.95, "Upper", color=color.red, linestyle=hline.style_dashed)
hline(0.05, "Lower", color=color.green, linestyle=hline.style_dashed)