norminv
// 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 An inverse normal distribution is a way to work backwards 
//	 from a known probability to find an x-value. It is an informal term and 
//   doesn't refer to a particular probability distribution. Returns the 
// 	 value of the inverse normal distribution function for a specified value, 
// 	 mean, and standard deviation.
// reference:
// 		https://github.com/teqniqly/norm-inv/blob/master/src/norm-inv.ts
//		https://support.microsoft.com/en-us/office/norminv-function-87981ab8-2de0-4cb0-b1aa-e21d4cb879b8
library("norminv")

// Error Function
erf(float xin)=>
    float[] cof = array.new<float>(28, 0)

    array.set(cof, 0, -1.3026537197817094)
    array.set(cof, 1, 6.4196979235649026e-1)
    array.set(cof, 2, 1.9476473204185836e-2)
    array.set(cof, 3, -9.561514786808631e-3)
    array.set(cof, 4, -9.46595344482036e-4)
    array.set(cof, 5, 3.66839497852761e-4)
    array.set(cof, 6, 4.2523324806907e-5)
    array.set(cof, 7, -2.0278578112534e-5)
    array.set(cof, 8, -1.624290004647e-6)
    array.set(cof, 9, 1.303655835580e-6)
    array.set(cof, 10, 1.5626441722e-8)
    array.set(cof, 11, -8.5238095915e-8)
    array.set(cof, 12, 6.529054439e-9)
    array.set(cof, 13, 5.059343495e-9)
    array.set(cof, 14, -9.91364156e-10)
    array.set(cof, 15, -2.27365122e-10)
    array.set(cof, 16, 9.6467911e-11)
    array.set(cof, 17, 2.394038e-12)
    array.set(cof, 18, -6.886027e-12)
    array.set(cof, 19,  8.94487e-13)
    array.set(cof, 20, 3.13092e-13)
    array.set(cof, 21, -1.12708e-13)
    array.set(cof, 22, 3.81e-16)
    array.set(cof, 23, 7.106e-15)
    array.set(cof, 24, -1.523e-15)
    array.set(cof, 25, -9.4e-17)
    array.set(cof, 26, 1.21e-16)
    array.set(cof, 27, -2.8e-17)

    bool isneg = false
    float d = 0
    float dd = 0
    float x = xin

    if xin < 0
        x := -xin
        isneg := true

    float t = 2 / (2 + x)
    float ty = 4 * t - 2

    for j = array.size(cof) - 1 to 1  
        float tmp = d
        d := ty * d - dd + array.get(cof, j)
        dd := tmp

    float res = t * math.exp(-x * x + 0.5 * (array.get(cof, 0) + ty * d) - dd)
    out = isneg ? res - 1 : 1 - res
    out

// Complementary Error Function
erfc(float x)=>
    out = 1 - erf(x)
    out

// Inverse Complementary Error Function
erfcinv(float prob)=>
    var float TWOSQRTPI = 1.12837916709551257
    float out = 0
    if prob >= 2 
        out := -100
    if prob <= 0
        out := 100
    else 
        float pp = (prob < 1) ? prob : 2 - prob
        float t = math.sqrt(-2 * math.log(pp / 2))
        float x = -0.70711 * ((2.30753 + t * 0.27061) / (1 + t * (0.99229 + t * 0.04481)) - t)
        for j = 0 to 2 - 1 
            float err = erfc(x) - pp
            x += err / (TWOSQRTPI * math.exp(-x * x) - x * err)
        out := (prob < 1) ? x : -x
    out

// @function Returns the value of the inverse normal distribution function for a specified value, mean, and standard deviation.
// @param x float, The input to the normal distribution function.
// @param mean float, The mean (mu) of the normal distribution function
// @param stdev float, The standard deviation (sigma) of the normal distribution function.
// @returns float.
export norminv(float x, float mean, float stdev)=>
    var float SQRTTWO = 1.41421356237309505
    float out = 0
    out := -SQRTTWO * stdev * erfcinv(2 * x) + mean

// example usage
out = norminv(0.5, 65.5, 2.5)
plot(out)