// Licensed under the Apache License, Version 2.0
// © mihakralj
//@version=6
indicator("Binomial Distribution CDF (BINOMDIST)", "BINOMDIST", 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 Log of binomial coefficient C(n, i) = ln(n!) - ln(i!) - ln((n-i)!)
//@param n Number of trials
//@param i Number of successes
//@returns ln(C(n, i))
lnBinom(simple int n, int i) =>
    lnGamma(float(n + 1)) - lnGamma(float(i + 1)) - lnGamma(float(n - i + 1))

//@function Binomial Distribution CDF: P(X <= k) = sum_{i=0}^{k} C(n,i) * p^i * (1-p)^(n-i)
//@param p Probability of success per trial (0 <= p <= 1)
//@param n Number of trials
//@param k Threshold (compute P(X <= k))
//@returns CDF value in [0,1]
//@optimized Log-space summation avoids factorial overflow for large n
binomCdf(series float p, simple int n, simple int k) =>
    if p <= 0.0
        k >= 0 ? 1.0 : 0.0
    else if p >= 1.0
        k >= n ? 1.0 : 0.0
    else
        float lnP = math.log(p)
        float lnQ = math.log(1.0 - p)
        float cdf = 0.0
        int kk = math.min(k, n)
        for i = 0 to kk
            float lnTerm = lnBinom(n, i) + i * lnP + (n - i) * lnQ
            cdf += math.exp(lnTerm)
        math.min(cdf, 1.0)

//@function Computes Binomial Distribution CDF for a normalized price series
//@param source Series to transform
//@param period Lookback period for min-max normalization to [0,1] as probability p
//@param trials Number of Bernoulli trials (n)
//@param threshold Success threshold (k) — compute P(X <= k)
//@returns Binomial CDF value in [0,1]
binomdist(series float source, simple int period, simple int trials, simple int threshold) =>
    if period <= 0
        runtime.error("Period must be greater than 0")
    if trials <= 0
        runtime.error("Trials must be greater than 0")
    if threshold < 0
        runtime.error("Threshold must be non-negative")

    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 p = range > 0.0 ? (source - minVal) / range : 0.5

    binomCdf(p, trials, threshold)

// ---------- 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_trials = input.int(20, "Trials (n)", minval=1, maxval=1000, tooltip="Number of Bernoulli trials")
i_threshold = input.int(10, "Threshold (k)", minval=0, maxval=1000, tooltip="Compute P(X <= k)")

// Calculation
result = binomdist(i_source, i_period, i_trials, i_threshold)

// Plot
plot(result, "BINOMDIST", 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)