Oceń określony przedział rozkładu normalnego

Wiem, że brakuje nieco łatwej w obsłudze formuły dla CDF normalnej dystrybucji, ze względu na skomplikowaną funkcję błędu.

Zastanawiam się jednak, czy istnieje fajna formuła dla . Albo jakie może być przybliżenie tego najnowszego stanu techniki.$N(c_{-} \leq x < c_{+}| \mu, \sigma^2)$

To zależy dokładnie od tego, czego szukasz . Poniżej znajduje się kilka krótkich szczegółów i odniesień.

Znaczna część literatury na temat aproksymacji koncentruje się wokół funkcji

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

dla . Wynika to z faktu, że podana funkcja może zostać rozłożona na zwykłą różnicę powyższej funkcji (ewentualnie dostosowaną przez stałą). Do tej funkcji odwołuje się wiele nazw, w tym „górny ogon rozkładu normalnego”, „prawa całka normalna normalna” i „ funkcja Q Gaussa ”, żeby wymienić tylko kilka. Zobaczysz także przybliżenia do współczynnika Millsa , który wynosi R ( x ) = Q ( x $x > 0$$Q$ , gdzieφ(x)=(2π)-1/2e-x2/2jest Gaussa PDF.

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

Oto kilka referencji do różnych celów, którymi możesz być zainteresowany.

Obliczeniowe

Rzeczywistym standardem obliczania funkcji lub powiązanej uzupełniającej funkcji błędu jest$Q$

WJ Cody, Rational Chebyshev Approximations for the Error Function , Math. Komp. , 1969, s. 631--637.

Każda (szanująca się) implementacja korzysta z tego dokumentu. (MATLAB, R itp.)

„Proste” przybliżenia

Abramowitz i Stegun mają jeden oparty na wielomianowym rozszerzeniu transformacji wejścia. Niektóre osoby używają go jako przybliżenia „o wysokiej precyzji”. Nie podoba mi się to do tego celu, ponieważ źle się zachowuje w pobliżu zera. Na przykład, ich zbliżenie to nie wytworzeniem Q ( 0 ) = 1 / 2 , które zdaniem jest duży nie-nie. Czasami zdarzają się z tego powodu złe rzeczy .$\hat{Q}(0) = 1/2$

Borjesson i Sundberg podają proste przybliżenie, które działa całkiem dobrze w większości aplikacji, w których wymaga się tylko kilku cyfr precyzji. Absolutny błąd względny nie gorsza niż 1%, co jest bardzo dobre biorąc pod uwagę jego prostota. Podstawowym aproksymacji P ( x ) = 1 i ich korzystne wybory stałe są=0,339ib=5,51. To odniesienie jest

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

PO Borjesson i CE Sundberg. Proste przybliżenia funkcji błędu Q (x) dla aplikacji komunikacyjnych . IEEE Trans. Commun , COM-27 (3): 639–643, marzec 1979 r.

Oto wykres jego bezwzględnego błędu względnego.

Literatura elektrotechniczna jest wypełniona różnymi tego rodzaju przybliżeniami i wydaje się, że nadmiernie się nimi interesuje. Wiele z nich jest biednych lub rozwija się w bardzo dziwne i skomplikowane wyrażenia.

Możesz także spojrzeć na

W. Bryc. Jednolite przybliżenie do prawej całki normalnej . Applied Mathematics and Computation , 127 (2-3): 365–374, kwiecień 2002.

Ciągła frakcja Laplace'a

$x > 0$

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

$1/(x+1/(x+2/(x+3/(x+\cdots))))$$x$$x = 0$

$Q(x)$

$$\frac{x}{x^2 + 1} < R(x) < \frac{1}{x} \>,$$ $Q$ $$\frac{x}{x^2 + 1} \varphi(x) < Q(x) < \frac{1}{x} \varphi(x) .$$ $x^{-2}$

$Q(x) \sim \varphi(x)/x$$a \in [0,1]$$x \to \infty$$b$

$Q$

$x$

CI. C. Lee. Na Laplace'u kontynuował ułamek normalny . Ann. Inst. Statystyk. Matematyka , 44 (1): 107–120, marzec 1992 r.

---

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

Mam nadzieję, że to zacznie. Jeśli masz bardziej konkretne zainteresowania, mogę cię gdzieś wskazać.

Przypuszczam, że spóźniłem się z bohaterem, ale chciałem skomentować post kardynała, a ten komentarz stał się zbyt duży dla zamierzonego pudełka.

$x >0$$x$

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

Istnieją w rzeczywistości alternatywne sposoby obliczania (uzupełniającej) funkcji błędu oprócz stosowania aproksymacji Czebyszewa. Ponieważ zastosowanie aproksymacji Czebyszewa wymaga przechowywania nie kilku współczynników, metody te mogą mieć przewagę, jeśli struktury tablic są nieco kosztowne w środowisku komputerowym (można wstawić współczynniki, ale wynikowy kod prawdopodobnie wyglądałby jak barok bałagan).

$|x|$, Abramowitz i Stegun dają ładnie zachowane serie (przynajmniej lepiej zachowane niż zwykła seria Maclaurin):

$$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}$$ (na podstawie wzoru 7.1.6 )

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

---

$|x|$

Lentz , Thompson i Barnett opracowali algorytm do liczbowej oceny ciągłej frakcji jako produktu nieskończonego, który jest bardziej wydajny niż zwykłe podejście do obliczania ciągłej frakcji „wstecz”. Zamiast wyświetlać ogólny algorytm, pokażę, w jaki sposób specjalizuje się on w obliczaniu współczynnika Millsa:

$\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}$

CF jest użyteczny tam, gdzie poprzednio wspomniana seria zaczyna się powoli zbiegać; będziesz musiał eksperymentować z określeniem odpowiedniego „punktu przerwania”, aby przełączyć się z serii na CF w środowisku komputerowym. Istnieje również alternatywa użycia serii asymptotycznej zamiast Laplaciana CF, ale z mojego doświadczenia wynika, że ​​Laplacian CF jest wystarczający do większości zastosowań.

---

Wreszcie, jeśli nie trzeba obliczyć funkcję (komplementarne) Błąd bardzo dokładnie (czyli tylko kilka cyfr znaczących), istnieją kompaktowe przybliżenia spowodowane Serge Winitzki. Oto jeden z nich:

$$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$

(Ta odpowiedź pierwotnie pojawiła się w odpowiedzi na podobne pytanie, następnie zamknięta jako duplikat. OP chciał tylko „implementacji” całki Gaussa, niekoniecznie „najnowszego stanu techniki”. W jego komentarzach stało się jasne, że stosunkowo prosta , preferowane byłoby krótkie wdrożenie).

---

Jak wskazują komentarze, musisz zintegrować plik PDF . Istnieje wiele sposobów wykonania całki. Dawno temu, kiedy obliczenia były powolne i kosztowne, David Hill opracował przybliżenie za pomocą prostej arytmetyki (funkcje wymierne i potęgowanie). Ma podwójną dokładność dla typowych argumentów (pomiędzy$-8.5$ i $+8.5$, w przybliżeniu). W 1973 roku opublikował wersję Fortran w Applied Statistics o nazwie ALNORM.F. Przez lata przenosiłem to do różnych środowisk, które nie miały całki normalnej (gaussowskiej) lub które miały podejrzane (np. Excel).

Wersja MatLab (z odpowiednimi atrybutami) jest dostępna na stronie http://people.sc.fsu.edu/~jburkardt/m_src/asa005/alnorm.m . Całkowicie nieudokumentowana wersja oryginalnego kodu Fortran pojawia się na stronie „Codeers Code Search” (sic).

Wiele lat temu przeniosłem to do AWK. Ta wersja może być bardziej odpowiednia dla współczesnego programisty do przenoszenia ze względu na jej składnię podobną do C (zamiast Fortran) i kilka dodatkowych komentarzy, które wstawiłem podczas jej opracowywania i testowania, ponieważ musiałem zwiększyć jej dokładność. Pojawia się poniżej.

Dla osób bez dużego doświadczenia w przenoszeniu kodu naukowego / matematycznego / statystycznego, kilka rad : jeden pojedynczy błąd typograficzny może spowodować poważne błędy, które mogą nie być łatwo wykryte. (Zaufaj mi, zrobiłem ich wiele.) Zawsze, zawsze wykonuj dokładny i wyczerpujący test. Ponieważ normalna funkcja całkowa / gaussowska całka / błąd jest dostępna w tak wielu tabelach i tak wielu programach, łatwo i szybko zestawia się w tabele ogromną liczbę wartości przenoszonej funkcji i systematycznie porównuje (tj. Z komputerem, nie na oko) wartości do poprawienia. Taki test można zobaczyć na początku mojego kodu: tworzy on tabelę wartości w -8,5: 8,5 (o 0,1), którą można przesłać (poprzez STDOUT) do innego programu w celu systematycznego sprawdzania.

Innym podejściem do testowania - dla tych, którzy mają wystarczającą wiedzę na temat analizy numerycznej, aby wiedzieć, jak oszacować oczekiwane błędy - byłoby numeryczne różnicowanie wartości i porównanie ich z plikiem PDF (który można łatwo obliczyć).

Nawiasem mówiąc: ten kod dotyczy tylko przypadku $0$i odchylenie standardowe jednostki („sigma”). Ale to wszystko, czego potrzebujesz: od integracji$-\infty$ do $x$ kiedy średnia to $\mu$ a SD jest $\sigma$, po prostu oblicz $z = (x-\mu)/\sigma$i zastosuj się alnormdo tego.

Edytować

Przetestowałem do portu alnormdo Mathematica, który oblicza wartości arbitralnej precyzji. Aby porównać wyniki, oto wykres naturalnego logu stosunków wartości górnych ogonów$1 - \Phi(z)$ z $z\ge 1$. (Dodatni błąd względny oznacza, że alnormjest zbyt duży).

Wartości są zawsze dokładne$4 \times 10^{-11}$ w odniesieniu do znikomo małych prawdopodobieństw ogona . Możesz zobaczyć, gdzie obliczenia przełączają się na asymptotyczną formułę (w$z=16$) i jest oczywiste, że ta formuła staje się niezwykle dokładna jako $z$wzrasta. Fabuła kończy się o godz$z=\sqrt(2\times 708) \approx 37.6$ ponieważ tutaj właśnie potęgowanie z podwójną precyzją zaczyna się przepełniać.

Na przykład alnorm[-6.0]zwraca$9.865\ 876\ 450\ 315E-10$ podczas gdy prawdziwa wartość, równa $\frac{1}{2}\text{erfc}(3\sqrt{2})$, jest w przybliżeniu $9.865\ 876\ 450\ 377E-10$, najpierw różniący się dwunastą cyfrą dziesiętną.

NB W ramach tej edycji, zmieniłem UPPER_TAIL_IS_ZEROod 15.celu 16.w kodzie: to sprawia, że wynik odrobinę bardziej dokładne dla$Z$ pomiędzy $15$ i $16$. (Koniec edycji.)

#----------------------------------------------------------------------#
#   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 ###