XProfan

Deutsch

p.specht

Die Umkehrfunktion für die einseitige Summenverteilung der Gauss'schen Glockenfunktion weist in Tabellen eine Genauigkeit von 6-7 Stellen auf. Es geht noch genauer, solange die jeweilige Maschinenpräzision mitspielt: Die nachstehende Variante ist durch Aufteilung auf zwei unterschiedliche Reihenentwicklungen für die Bereiche um die Mitte und an den Flügeln der Glockenkurve garantiert auf 9 Stellen genau, also im ppm-(parts per million)-Bereich. Und trotzdem eine halbwegs flotte Demo!

Kompilieren Markieren Separieren

WindowTitle "Inverse Normalverteilung"
' Fortran Source: ren-raw chen, rutgers business school: Normal Distribution Inversed
' Transposed to XProfan-11.2a (CL) 2014-10 by P. Specht, Wien
' No Warranty whatsoever! Keine wie auch immer geartete Gewähr!
WindowStyle 24:set("decimals",17):declare p$:CLS
doitagain:
Print " Gewünschte Eintrittswahrscheinlichkeit [%]?: ";:input p$:case p$="":end
print InvNormDist(val(p$)/100);" * Sigma":if %csrlin>40:waitinput:cls:endif
goto "doitagain"

proc InvNormDist :parameters p!

    if p!<0:Print "*** Error: Negative Wahrscheinlichkeiten?***":return -999999:beep:endif

        if p!>0.99999999999999994:Print "*** Abs.Sicherheit nur bei Abdeckung bis +Inf!***":return 10^53:endif

            if p!<val("1E-53"):Print "*** Sicher kein Eintritt? Nur bei -Inf!***":return val("-1E53"):endif

                ' Source: ren-raw chen, rutgers business school: ' normal distribution inversed
                ' (translated from https://home.online.no/~pjacklam/notes/invnorm by john herrero)
                declare p_low!,p_high!, z!,q!,r!
                declare a1!,a2!,a3!,a4!,a5!,a6!, b1!,b2!,b3!,b4!,b5!
                declare c1!,c2!,c3!,c4!,c5!,c6!, d1!,d2!,d3!,d4!
                a1!=-39.6968302866538:a2!=220.946098424521:a3!=-275.928510446969:a4!=138.357751867269
                a5!=-30.6647980661472:a6!=2.50662827745924:b1!=-54.4760987982241:b2!=161.585836858041
                b3!=-155.698979859887:b4!=66.8013118877197:b5!=-13.2806815528857:c1!=-0.00778489400243029
                c2!=-0.322396458041136:c3!=-2.40075827716184:c4!=-2.54973253934373
                c5!=4.37466414146497:c6!=2.93816398269878:d1!=0.00778469570904146
                d2!=0.32246712907004:d3!=2.445134137143:d4!=3.75440866190742:p_low!=0.02425:p_high!=1-p_low!
                case p! <p_low!: goto "g201"
                case p!>=p_low!: goto "g301"
                g201:
                q!=sqrt(-2*ln(p!))
                z!=(((((c1!*q!+c2!)*q!+c3!)*q!+c4!)*q!+c5!)*q!+c6!)/ \
                ((((d1!*q!+d2!)*q!+d3!)*q!+d4!)*q!+1)
                goto "g204"
                g301:
                case (p!>=p_low!) and (p!<p_high!): goto "g202"
                case p!>p_high!: goto "g302"
                g202:
                q!=p!-0.5
                r!=q!*q!
                z!=(((((a1!*r!+a2!)*r!+a3!)*r!+a4!)*r!+a5!)*r!+a6!)*q! / \
                (((((b1!*r!+b2!)*r!+b3!)*r!+b4!)*r!+b5!)*r!+1)
                goto "g204"
                g302:
                case (p!>p_high!) and (p!<1): goto "g203"
                g203:
                q!=sqrt(-2*ln(1-p!))
                z!=-(((((c1!*q!+c2!)*q!+c3!)*q!+c4!)*q!+c5!)*q!+c6!) / \
                ((((d1!*q!+d2!)*q!+d3!)*q!+d4!)*q!+1)
                g204:
                return z!

            endproc