XProfan

Windowtitle "Approximation der Inversen Funktion des Normalverteilungsintegrals "
' (D) Demo by P. Specht 07-2013, ohne jegliche Gewähr!
Windowstyle 1048:Cls:font 2
Declare p!

Repeat

    Print "\n Gewünschte symmetrische StandardNormal-\n verteilte Wahrscheinlichkeit [%]:?",
    input p!
    case %csrlin>24:cls
    p!=abs(p!)
    case p!>100:p!=100
    p!=p!/100
    Print " Ergebnis: ";format$("%g", SNormInv(p!) )
    print "-------------------------------------------------"

until 0

End

proc SNormInv :parameters p!

    '  Acklam's InverseCDF 'Inverse Normalverteilung
    '  Adapted for XProfan 11.2a from Peter Acklam's
    '  "An algorithm for computing the inverse normal cumulative distribution function"
    '  (https://home.online.no/~pjacklam/notes/invnorm/)
    '  by P. Specht (July 26th, 2013, adding conditional definition of constants)
    '  Coefficients for rational approximations (Padé?)
    declare q!,r!
    'Define break-points:
    var P_LOW!= 0.02425
    var P_HIGH!= 1 - P_LOW!

    If (p!>0) And (p!<P_LOW!)

        'Rational approximation for lower region.
        q!=Sqrt(-2*Lg(p!))
        var C1! = -0.0077848940024302926
        var C2! = -0.32239645804113648
        var C3! = -2.4007582771618381
        var C4! = -2.5497325393437338
        var C5! = 4.3746641414649678
        var C6! = 2.9381639826987831
        var D1! = 0.0077846957090414622
        var D2! = 0.32246712907003983
        var D3! = 2.445134137142996
        var D4! = 3.7544086619074162
        Return (((((C1!*q!+C2!)*q!+C3!)*q!+C4!)*q!+C5!)*q!+C6!)/ \
        ((((D1!*q!+D2!)*q!+D3!)*q!+D4!)*q!+1)

    ElseIf (p!>=P_LOW!) And (p!<=P_HIGH!)

        'Rational approximation for central region
        q!=p!-0.5
        r!=q!*q!
        var A1! = -39.696830286653757
        var A2! = 220.9460984245205
        var A3! = -275.92851044696869
        var A4! = 138.357751867269
        var A5! = -30.66479806614716
        var A6! = 2.5066282774592392
        var B1! = -54.476098798224058
        var B2! = 161.58583685804089
        var B3! = -155.69897985988661
        var B4! = 66.80131188771972
        var B5! = -13.280681552885721
        Return (((((A1!*r!+A2!)*r!+A3!)*r!+A4!)*r!+A5!)*r!+A6!)*q!/ \
        (((((B1!*r!+B2!)*r!+B3!)*r!+B4!)*r!+B5!)*r!+1)

    ElseIf (p!>P_HIGH!) And (p!<1)

        'Rational approximation for upper region
        q!=Sqrt(-2*Lg(1-p!))
        ' Definition wurde in den IF-Teil verlegt, daher hier doppelt definiert (sh. oben)
        var C1! = -0.0077848940024302926
        var C2! = -0.32239645804113648
        var C3! = -2.4007582771618381
        var C4! = -2.5497325393437338
        var C5! = 4.3746641414649678
        var C6! = 2.9381639826987831
        var D1! = 0.0077846957090414622
        var D2! = 0.32246712907003983
        var D3! = 2.445134137142996
        var D4! = 3.7544086619074162
        Return -1*(((((C1!*q!+C2!)*q!+C3!)*q!+C4!)*q!+C5!)*q!+C6!)/ \
        ((((D1!*q!+D2!)*q!+D3!)*q!+D4!)*q!+1)

    Else

        Print " ERROR: Argument Out Of Range! "
        beep
        sound 2000,200
        return val("-1e-37")

    EndIf

Endproc