Bestimmte Intervalle der Normalverteilung auswerten

Ich weiß, dass eine einfach zu handhabende Formel für die CDF einer Normalverteilung aufgrund der darin enthaltenen komplizierten Fehlerfunktion etwas fehlt.

Ich frage mich jedoch, ob es eine schöne Formel für $N(c_{-} \leq x < c_{+}| \mu, \sigma^2)$ . Oder was die "State-of-the-Art" -Näherung für dieses Problem sein könnte.

Es kommt genau darauf an, wonach Sie suchen . Nachfolgend finden Sie einige kurze Details und Referenzen.

Ein Großteil der Literatur für Approximationen konzentriert sich auf die Funktion

$$Q(x) = \int_x^\infty \frac{1}{\sqrt{2\pi}} e^{-\frac{u^2}{2}} \, \mathrm{d}u$$

für . Dies liegt daran, dass die von Ihnen bereitgestellte Funktion als einfache Differenz der obigen Funktion zerlegt werden kann (möglicherweise durch eine Konstante angepasst). Auf diese Funktion wird mit vielen Namen Bezug genommen, einschließlich "oberes Ende der Normalverteilung", "rechtes normales Integral" und "Gaußsche ", um nur einige zu nennen. Außerdem sehen Sie Annäherungen an das Mills-Verhältnis : wobei ist das Gaußsche pdf.$x > 0$$Q$

$$R(x) = \frac{Q(x)}{\varphi(x)}$$ $\varphi(x) = (2\pi)^{-1/2} e^{-x^2 / 2}$

Hier liste ich einige Referenzen für verschiedene Zwecke auf, die Sie interessieren könnten.

Computational

Der De-facto-Standard zur Berechnung der Funktion oder der damit verbundenen komplementären Fehlerfunktion ist$Q$

WJ Cody, Rational Chebyshev Approximationen für die Fehlerfunktion , Math. Comp. 1969, S. 631-637.

Jede (selbst-respektierende) Implementierung verwendet dieses Papier. (MATLAB, R usw.)

"Einfache" Annäherungen

Abramowitz und Stegun haben eine polynomielle Erweiterung einer Transformation der Eingabe zugrunde gelegt. Einige Leute benutzen es als "hochpräzise" Näherung. Aus diesem Grund mag ich es nicht, da es sich bei Null schlecht verhält. Zum Beispiel ihre Annäherung ist nicht nachgeben , was meiner Meinung nach ein großes No-No. Manchmal passieren daraus schlimme Dinge .$\hat{Q}(0) = 1/2$

Borjesson und Sundberg geben eine einfache Näherung an, die sich für die meisten Anwendungen eignet, bei denen nur wenige Stellen Genauigkeit erforderlich sind. Der absolute relative Fehler ist niemals schlechter als 1%, was angesichts seiner Einfachheit recht gut ist. Der Grund Näherung ist Q ( x ) = 1 und deren bevorzugte Wahl der Konstanten sinda=0,339undb=5,51. Diese Referenz ist

$$\hat{Q}(x) = \frac{1}{(1-a) x + a \sqrt{x^2 + b}} \varphi(x)$$ $a = 0.339$$b = 5.51$

PO Borjesson und CE Sundberg. Einfache Näherungen der Fehlerfunktion Q (x) für Kommunikationsanwendungen . IEEE Trans. Kommun. , COM-27 (3): 639–643, März 1979.

Hier ist eine Darstellung des absoluten relativen Fehlers.

Die elektrotechnische Literatur steckt voller solcher Näherungen und scheint sich übermäßig intensiv mit ihnen zu beschäftigen. Viele von ihnen sind arm oder zeigen sehr merkwürdige und verschlungene Ausdrücke.

Sie könnten sich auch anschauen

W. Bryc. Eine gleichmäßige Annäherung an das rechte normale Integral . Applied Mathematics and Computation , 127 (2-3): 365–374, April 2002.

Laplace's fortgesetzte Fraktion

Laplace hat eine schöne fortgesetzte Fraktion, die aufeinanderfolgende obere und untere Schranken für jeden Wert von ergibt . Es ist, ausgedrückt in Mills 'Ratio,$x > 0$

$$R(x) = \frac{1}{x+}\frac{1}{x+}\frac{2}{x+}\frac{3}{x+}\cdots ,$$

wobei die von mir verwendete Notation für einen fortgesetzten Bruch ziemlich normal ist , dh . Dieser Ausdruck konvergiert jedoch für kleines x nicht sehr schnell und divergiert bei x = 0 .$1/(x+1/(x+2/(x+3/(x+\cdots))))$$x$$x = 0$

Diese fortgesetzte Fraktion liefert tatsächlich viele der "einfachen" Schranken für , die Mitte bis Ende des 20. Jahrhunderts "wiederentdeckt" wurden. Es ist leicht zu erkennen, dass für einen fortgesetzten Bruch in "Standard" -Form (dh bestehend aus positiven ganzzahligen Koeffizienten) das Abschneiden des Bruches bei ungeraden (geraden) Termen eine obere (untere) Grenze ergibt.$Q(x)$

Laplace sagt uns daher sofort, dass Beide sind Grenzendie in der Mitte der 1900er Jahre „neu entdeckt“ wurden. In Bezug auf die Q-Funktion entspricht dies

$$\frac{x}{x^2 + 1} < R(x) < \frac{1}{x} \>,$$ $Q$ Ein alternativer Beweis dafür durch einfache Teilintegration findet sich in S. Resnick,Adventures in Stochastic Processes, Birkhauser, 1992, in Kapitel 6 (Brownsche Bewegung). Der absolute relative Fehler dieser Grenzen ist nicht schlechter alsx-2, wie indieser verwandten Antwort gezeigt. $$\frac{x}{x^2 + 1} \varphi(x) < Q(x) < \frac{1}{x} \varphi(x) .$$ $x^{-2}$

Beachten Sie insbesondere, dass die obigen Ungleichungen unmittelbar bedeuten, dass . Diese Tatsache kann auch mit der Regel von L'Hopital festgestellt werden. Dies erklärt auch die Wahl der funktionalen Form der Borjesson-Sundberg-Näherung. Jede Wahl eines ∈ [ 0 , 1 ] hält die asymptotische Äquivalenz als x → ∞ . Der Parameter b dient als "Durchgangskorrektur" nahe Null.$Q(x) \sim \varphi(x)/x$$a \in [0,1]$$x \to \infty$$b$

Hier ist eine Darstellung der und der beiden Laplace-Grenzen.$Q$

$x$

CI. C. Lee. Auf Laplace wird der Bruch für das normale Integral fortgesetzt . Ann. Inst. Statist. Mathematik. 44 (1), 107–120 (März 1992).

---

$Q(x)$$x$$x > 3$

Hoffentlich können Sie damit anfangen. Wenn Sie ein spezifischeres Interesse haben, kann ich Sie möglicherweise auf etwas hinweisen.

Ich denke, ich bin zu spät, der Held, aber ich wollte Kardinals Beitrag kommentieren, und dieser Kommentar wurde zu groß für die vorgesehene Box.

$x >0$$x$

$\mathrm{erf}(x)$$R(x)$

Abgesehen von der Verwendung von Chebyshev-Näherungen gibt es in der Tat alternative Möglichkeiten zur Berechnung der (komplementären) Fehlerfunktion. Da für die Verwendung einer Chebyshev-Näherung nicht wenige Koeffizienten gespeichert werden müssen, haben diese Methoden möglicherweise einen Vorteil, wenn Array-Strukturen in Ihrer Computerumgebung etwas kostspielig sind (Sie könnten die Koeffizienten einbetten, aber der resultierende Code würde wahrscheinlich wie ein Barock aussehen Chaos).

$|x|$

$$R(x)=\sqrt{\frac{\pi}{2}}\exp\left(\frac{x^2}{2}\right)-x\sum_{j=0}^\infty\frac{2^j j!}{(2j+1)!}x^{2j}$$

$x^{2j}$$c_j=\frac{2^j j!}{(2j+1)!}$$c_0=1$$c_{j+1}=\frac{c_j}{2j+3}$

---

$|x|$

Lentz , Thompson und Barnett haben einen Algorithmus zur numerischen Bewertung einer fortgesetzten Fraktion als unendliches Produkt abgeleitet, der effizienter ist als der übliche Ansatz, eine fortgesetzte Fraktion "rückwärts" zu berechnen. Anstatt den allgemeinen Algorithmus anzuzeigen, werde ich zeigen, wie er sich auf die Berechnung des Mills-Verhältnisses spezialisiert:

$\displaystyle Y_0=x,\,C_0=Y_0,\,D_0=0$
$\text{repeat for }j=1,2,\dots$

$$D_j=\frac1{x+jD_{j-1}}$$ $$C_j=x+\frac{j}{C_{j-1}}$$ $$H_j=C_j D_j$$ $$Y_j=H_j Y_{j-1}$$ $\text{until }|H_j-1| < \text{tol}$
$\displaystyle R(x)=\frac1{Y_j}$

$\text{tol}$

Die CF ist nützlich, wenn die zuvor erwähnte Serie langsam konvergiert. Sie müssen experimentieren, um den geeigneten "Haltepunkt" für den Wechsel von der Serie zur CF in Ihrer Computerumgebung zu bestimmen. Es gibt auch die Alternative, anstelle des Laplace-CF eine asymptotische Serie zu verwenden, aber meiner Erfahrung nach ist der Laplace-CF für die meisten Anwendungen gut genug.

---

Wenn Sie die (komplementäre) Fehlerfunktion nicht sehr genau berechnen müssen (dh nur auf wenige signifikante Stellen genau), gibt es kompakte Näherungen, die auf Serge Winitzki zurückzuführen sind. Hier ist einer von ihnen:

$$R(x)\approx \frac{\sqrt{2\pi}+x(\pi-2)}{2+x\sqrt{2\pi}+x^2(\pi-2)}$$

$1.84\times 10^{-2}$$x$

(Diese Antwort erschien ursprünglich als Antwort auf eine ähnliche Frage, die später als Duplikat geschlossen wurde. Das OP wollte nur "eine" Implementierung des Gaußschen Integrals, nicht notwendigerweise "Stand der Technik". In seinen Kommentaren wurde deutlich, dass eine relativ einfache wäre eine kurze Implementierung vorzuziehen.)

---

$-8.5$$+8.5$

Eine MatLab-Version (mit entsprechenden Attributen) ist unter http://people.sc.fsu.edu/~jburkardt/m_src/asa005/alnorm.m verfügbar . Eine vollständig undokumentierte Version des ursprünglichen Fortran-Codes wird auf einer "Koders Code Search" -Site (sic) angezeigt .

Vor vielen Jahren habe ich das auf AWK portiert. Diese Version ist für den modernen Entwickler aufgrund der C-ähnlichen Syntax (anstelle von Fortran) und einiger zusätzlicher Kommentare, die ich beim Entwickeln und Testen eingefügt habe, möglicherweise günstiger, da ich die Genauigkeit verbessern musste. Es erscheint unten.

Für diejenigen, die nicht viel Erfahrung mit dem Portieren von wissenschaftlichem / mathematischem / Statistik-Code haben, einige Ratschläge : Ein einziger Tippfehler kann schwerwiegende Fehler verursachen, die möglicherweise nicht leicht zu erkennen sind. (Vertrauen Sie mir, ich habe viele davon gemacht.) Machen Sie immer einen sorgfältigen und umfassenden Test. Da die normale Integral- / Gaußsche Integral- / Fehlerfunktion in so vielen Tabellen und in so vielen Programmen verfügbar ist, können Sie einfach und schnell eine große Anzahl von Werten Ihrer portierten Funktion tabellieren und systematisch vergleichen (dh mit dem Computer, nicht mit dem Auge). die zu korrigierenden Werte. Sie können einen solchen Test am Anfang meines Codes sehen: Er erzeugt eine Wertetabelle in -8,5: 8,5 (mal 0,1), die (über STDOUT) an ein anderes Programm zur systematischen Überprüfung weitergeleitet werden kann.

Ein anderer Testansatz - für diejenigen mit genügend numerischem Hintergrund, um die zu erwartenden Fehler abzuschätzen - wäre die numerische Unterscheidung der Werte und der Vergleich mit dem PDF (das leicht berechnet werden kann).

$0$$-\infty$$x$$\mu$$\sigma$$z = (x-\mu)/\sigma$alnorm

Bearbeiten

alnorm$1 - \Phi(z)$$z\ge 1$alnorm

$4 \times 10^{-11}$ $z=16$$z$$z=\sqrt(2\times 708) \approx 37.6$

alnorm[-6.0]$9.865\ 876\ 450\ 315E-10$$\frac{1}{2}\text{erfc}(3\sqrt{2})$$9.865\ 876\ 450\ 377E-10$

UPPER_TAIL_IS_ZERO15.16.$Z$$15$$16$

#----------------------------------------------------------------------#
#   ALNORM.AWK
#   Compute values of the cumulative normal probability function.
#   From G. Dallal's STAT-SAK (Fortran code).
#   Additional precision using asymptotic expression added 7/8/92.
#----------------------------------------------------------------------#
BEGIN {
    for (i=-85; i<=85; i++) {
        x = i/10
        p = alnorm(x, 0)
        printf("%3.1f %12.10f\n", x, p)
    }
    exit
}
function alnorm(z,up,    y,aln,w) {
#
#    ALGORITHM AS 66 APPL. STATIST. (1973) VOL.22, NO.3:
#    Hill,  I.D.  (1973).  Algorithm AS 66.  The normal  integral.
#                          Appl. Statist.,22,424-427.
#
#    Evaluates the tail area of the standard normal curve from
#    z to infinity if up, or from -infinity to z if not up.
#
#    LOWER_TAIL_IS_ONE, UPPER_TAIL_IS_ZERO, and EXP_MIN_ARG
#    must be set to suit this computer and compiler.

    LOWER_TAIL_IS_ONE = 8.5     # I.e., alnorm(8.5,0) = .999999999999+
    UPPER_TAIL_IS_ZERO = 16.0   # Changes to power series expression
    FORMULA_BREAK = 1.28        # Changes cont. fraction coefficients
    EXP_MIN_ARG = -708          # I.e., exp(-708) is essentially true 0

    if (z < 0.0) {
        up = !up
        z = -z
    }
    if ((z <= LOWER_TAIL_IS_ONE) || (up && z <= UPPER_TAIL_IS_ZERO)) {
        y = 0.5 * z * z
        if (z > FORMULA_BREAK) {
            if (-y > EXP_MIN_ARG) {
                aln = .398942280385 * exp(-y) / \
                  (z - 3.8052E-8 + 1.00000615302 / \
                  (z + 3.98064794E-4 + 1.98615381364 / \
                  (z - 0.151679116635 + 5.29330324926 / \
                  (z + 4.8385912808 - 15.1508972451 / \
                  (z + 0.742380924027 + 30.789933034 / \
                  (z + 3.99019417011))))))
            } else {
                aln = 0.0
            }
        } else {
            aln = 0.5 - z * (0.398942280444 - 0.399903438504 * y / \
              (y + 5.75885480458 - 29.8213557808 / \
              (y + 2.62433121679 + 48.6959930692 / \
              (y + 5.92885724438))))
        }
    } else {
        if (up) {   # 7/8/92
            # Uses asymptotic expansion for exp(-z*z/2)/alnorm(z)
            # Agrees with continued fraction to 11 s.f. when z >= 15
            # and coefficients through 706 are used.
            y = -0.5*z*z
            if (y > EXP_MIN_ARG) {
                w = -0.5/y  # 1/z^2
                aln = 0.3989422804014327*exp(y)/ \
                    (z*(1 + w*(1 + w*(-2 + w*(10 + w*(-74 + w*706))))))
                # Next coefficients would be -8162, 110410
            } else {
                aln = 0.0
            }
        } else {
            aln = 0.0
        }
    }
    return up ? aln : 1.0 - aln
}
### end of file ###