chi2Inv
// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © loxx

//@version=5
// @description Returns the inverse of the right-tailed or left-tailed probability of the chi-squared distribution.
// reference:
//      https://github.com/zetos/inv-chisquare-cdf

library("chi2Inv")

// import gamma log functions
import RicardoSantos/MathSpecialFunctionsGamma/1 as gl

//inverse regularized low gamma function
invRegLowGamma(p, a)=>
    float a1 = a - 1
    float EPS = 1e-8
    float gln = gl.GammaLn(a)
    float inverseRegLowGamma = 0
    float err = 0
    float t = 0
    float u = 0
    float pp = 0
    float lna1 = 0
    float afac = 0
    float out = 0

    if a > 1
        lna1 := math.log(a1)
        afac := math.exp(a1 * (lna1 - 1) - gln)
        pp := (p < 0.5) ? p : 1 - p
        t := math.sqrt(-2 * math.log(pp))
        inverseRegLowGamma := (2.30753 + t * 0.27061) / (1 + t * (0.99229 + t * 0.04481)) - t
        if (p < 0.5) 
            inverseRegLowGamma := -inverseRegLowGamma
        inverseRegLowGamma := math.max(1e-3, a * math.pow(1 - 1 / (9 * a) - inverseRegLowGamma / (3 * math.sqrt(a)), 3))
    else 
        t := 1 - a * (0.253 + a * 0.12)
        if p < t
            inverseRegLowGamma := math.pow(p / t, 1 / a)
        else 
            inverseRegLowGamma := 1 - math.log(1 - (p - t) / (1 - t))
    for j = 0 to 12 - 1
        if inverseRegLowGamma <= 0
            out := 0
        err := gl.GammaLowerRegularized(a, inverseRegLowGamma) - p
        if a > 1
            t := afac * math.exp(-(inverseRegLowGamma - a1) + a1 * (math.log(inverseRegLowGamma) - lna1))
        else 
            t := math.exp(-inverseRegLowGamma + a1 * math.log(inverseRegLowGamma) - gln)
        u := err / t
        t := u / (1 - 0.5 * math.min(1, u * ((a - 1) / inverseRegLowGamma - 1)))
        inverseRegLowGamma -= t
        if inverseRegLowGamma <= 0
            inverseRegLowGamma := 0.5 * (inverseRegLowGamma + t)
        if math.abs(t) < EPS * inverseRegLowGamma
            break

    out := inverseRegLowGamma

// @function Returns the inverse of the right-tailed or left-tailed probability of the chi-squared distribution.
// @param pin float, probability 
// @param deg float, degress of freedom.
// @param tail string, tail, either "r" for right or "l" for left
// @returns float.
export chi2Inv(float pin, float deg, string tail)=>
    float p = tail == "l" ? pin : 1 - pin
    out = 2 * invRegLowGamma(p, 0.5 * deg)  
    out

//example usage
right = chi2Inv(.93, 2, "r")
left = chi2Inv(.93, 2, "l")

plot(right, "right")
plot(left, "left")