Mengevaluasi interval pasti dari distribusi normal

Saya tahu bahwa formula CDF yang mudah ditangani untuk distribusi normal agak hilang, karena fungsi kesalahan yang rumit di dalamnya.

Namun, saya bertanya-tanya apakah ada formula yang bagus untuk $N(c_{-} \leq x < c_{+}| \mu, \sigma^2)$ . Atau apa kira-kira pendekatan "canggih" untuk masalah ini.

Tergantung pada apa yang Anda cari . Di bawah ini adalah beberapa perincian dan referensi singkat.

Banyak literatur untuk aproksimasi berpusat di sekitar fungsi

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

untuk $x > 0$ . Ini karena fungsi yang Anda berikan dapat didekomposisi sebagai perbedaan sederhana dari fungsi di atas (mungkin disesuaikan oleh konstanta). Fungsi ini disebut dengan banyak nama, termasuk "ekor atas dari distribusi normal", "integral normal kanan", dan "Gaussian $Q$ -function", untuk beberapa nama. Anda juga akan melihat perkiraan rasio Mills , yaitu di manaφ(x)=(2π)-1/2e-x2/2adalah Gaussian pdf.

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

Di sini saya daftar beberapa referensi untuk berbagai keperluan yang Anda mungkin tertarik.

Komputasi

Standar de-facto untuk menghitung fungsi- atau fungsi kesalahan komplementer terkait adalah$Q$

WJ Cody, Perkiraan Chebyshev Rasional untuk Fungsi Kesalahan , Matematika. Comp. , 1969, hlm. 631--637.

Setiap implementasi (menghargai diri sendiri) menggunakan tulisan ini. (MATLAB, R, dll.)

Perkiraan "Sederhana"

Abramowitz dan Stegun memiliki satu didasarkan pada ekspansi polinomial dari suatu transformasi input. Beberapa orang menggunakannya sebagai perkiraan "presisi tinggi". Saya tidak suka untuk tujuan itu karena berperilaku buruk di sekitar nol. Sebagai contoh, pendekatan mereka tidak tidak menghasilkan Q ( 0 ) = 1 / 2 , yang saya pikir besar tidak-tidak. Terkadang hal buruk terjadi karena ini.$\hat{Q}(0) = 1/2$

Borjesson dan Sundberg memberikan perkiraan sederhana yang bekerja cukup baik untuk sebagian besar aplikasi di mana seseorang hanya membutuhkan beberapa digit presisi. The kesalahan relatif absolut tidak pernah lebih buruk dari 1%, yang cukup baik mengingat kesederhanaan. Pendekatan dasar adalah Q ( x ) = 1 dan pilihan konstanta pilihannya adalaha=0,339danb=5,51. Referensi itu

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

PO Borjesson dan CE Sundberg. Perkiraan sederhana dari fungsi kesalahan Q (x) untuk aplikasi komunikasi . IEEE Trans. Komunal. , COM-27 (3): 639–643, Maret 1979.

Berikut adalah plot kesalahan relatif absolutnya.

Literatur teknik-listrik dipenuhi dengan berbagai perkiraan seperti itu dan tampaknya sangat tertarik pada mereka. Banyak dari mereka yang miskin atau berkembang menjadi ekspresi yang sangat aneh dan berbelit-belit.

Anda mungkin juga melihat

W. Bryc. Suatu pendekatan seragam ke integral normal yang tepat . Matematika dan Komputasi Terapan , 127 (2-3): 365-374, April 2002.

Fraksi lanjutan Laplace

Laplace memiliki fraksi lanjutan yang indah yang menghasilkan batas atas dan bawah berturut-turut untuk setiap nilai . Dalam hal rasio Mills,$x > 0$

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

di mana notasi yang saya gunakan cukup standar untuk fraksi lanjutan , yaitu, . Ekspresi ini tidak konvergen sangat cepat untuk x kecil , dan itu menyimpang di x = 0 .$1/(x+1/(x+2/(x+3/(x+\cdots))))$$x$$x = 0$

Fraksi lanjutan ini sebenarnya menghasilkan banyak batasan "sederhana" pada yang "ditemukan kembali" pada pertengahan hingga akhir 1900-an. Sangat mudah untuk melihat bahwa untuk fraksi lanjutan dalam bentuk "standar" (yaitu, terdiri dari koefisien bilangan bulat positif), memotong fraksi pada istilah ganjil (genap) memberikan batas atas (bawah).$Q(x)$

Karenanya, Laplace segera memberi tahu kita bahwa keduanya merupakan batas yang "ditemukan kembali" pada pertengahan 1900-an. Dalam halfungsi- Q , ini setara dengan

$$\frac{x}{x^2 + 1} < R(x) < \frac{1}{x} \>,$$ $Q$ Bukti alternatif dari ini menggunakan integrasi sederhana oleh bagian dapat ditemukan di S. Resnick, Adventures in Stochastic Processes , Birkhauser, 1992, dalam Bab 6 (gerak Brown). Kesalahan relatif absolut dari batas ini tidak lebih buruk dari x - 2 , seperti yang ditunjukkan dalamjawaban terkait ini. $$\frac{x}{x^2 + 1} \varphi(x) < Q(x) < \frac{1}{x} \varphi(x) .$$ $x^{-2}$

Perhatikan, khususnya, bahwa ketidaksetaraan di atas segera menyiratkan bahwa . Fakta ini dapat dibangun dengan menggunakan aturan L'Hopital juga. Ini juga membantu menjelaskan pilihan bentuk fungsional perkiraan Borjesson-Sundberg. Setiap pilihan yang ∈ [ 0 , 1 ] mempertahankan kesetaraan asimtotik sebagai x → ∞ . Parameter b berfungsi sebagai "koreksi kontinuitas" mendekati nol.$Q(x) \sim \varphi(x)/x$$a \in [0,1]$$x \to \infty$$b$

Berikut adalah plot fungsi- dan dua batas Laplace.$Q$

$x$

CI. C. Lee. On Laplace melanjutkan pecahan untuk integral normal . Ann. Inst. Statist. Matematika , 44 (1): 107–120, Maret 1992.

---

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

Semoga ini akan membantu Anda memulai. Jika Anda memiliki minat yang lebih spesifik, saya mungkin bisa mengarahkan Anda ke suatu tempat.

Saya kira saya sudah terlambat pahlawan, tetapi saya ingin mengomentari posting kardinal, dan komentar ini menjadi terlalu besar untuk kotak yang dimaksud.

$x >0$$x$

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

Sebenarnya ada cara alternatif untuk menghitung fungsi kesalahan (komplementer) selain dari menggunakan pendekatan Chebyshev. Karena penggunaan pendekatan Chebyshev memerlukan penyimpanan tidak sedikit koefisien, metode ini mungkin memiliki keunggulan jika struktur array sedikit mahal di lingkungan komputasi Anda (Anda dapat menyejajarkan koefisien, tetapi kode yang dihasilkan mungkin akan terlihat seperti barok kekacauan).

$|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 dan Barnett memperoleh algoritma untuk mengevaluasi secara numerik fraksi lanjutan sebagai produk tanpa batas, yang lebih efisien daripada pendekatan biasa menghitung fraksi lanjutan "mundur". Alih-alih menampilkan algoritma umum, saya akan menunjukkan bagaimana itu mengkhususkan diri pada perhitungan rasio Mills:

$\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 berguna di mana seri yang disebutkan sebelumnya mulai bertemu perlahan; Anda harus bereksperimen dengan menentukan "break point" yang tepat untuk beralih dari seri ke CF di lingkungan komputasi Anda. Ada juga alternatif menggunakan seri asimptotik alih-alih CF Laplacian, tetapi pengalaman saya adalah bahwa CF Laplacian cukup baik untuk sebagian besar aplikasi.

---

Akhirnya, jika Anda tidak perlu menghitung fungsi kesalahan (komplementer) dengan sangat akurat (yaitu hanya beberapa digit signifikan), ada perkiraan ringkas karena Serge Winitzki. Ini salah satunya:

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

Perkiraan ini memiliki kesalahan relatif maksimum $1.84\times 10^{-2}$ dan menjadi lebih akurat sebagai $x$ meningkat.

(Jawaban ini awalnya muncul sebagai jawaban atas pertanyaan serupa, kemudian ditutup sebagai duplikat. OP hanya menginginkan "sebuah" implementasi integral Gaussian, tidak harus "canggih". Dalam komentarnya, tampak jelas bahwa relatif sederhana , implementasi singkat akan lebih disukai.)

---

Seperti yang ditunjukkan oleh komentar, Anda perlu mengintegrasikan PDF . Ada banyak cara untuk melakukan integral. Dahulu kala, ketika perhitungannya lambat dan mahal, David Hill menyusun perkiraan menggunakan aritmatika sederhana (fungsi rasional dan eksponensial). Ini memiliki akurasi presisi ganda untuk argumen umum (antara$-8.5$ dan $+8.5$, kurang-lebih). Pada tahun 1973 ia menerbitkan versi Fortran dalam Statistik Terapan disebut ALNORM.F. Selama bertahun-tahun saya telah porting ini ke berbagai lingkungan yang tidak memiliki integral Normal (Gaussian) atau yang memiliki yang dicurigai (seperti Excel).

Versi MatLab (dengan atribusi yang sesuai) tersedia di http://people.sc.fsu.edu/~jburkardt/m_src/asa005/alnorm.m . Versi kode Fortran yang sepenuhnya tidak berdokumen muncul di situs "Koders Code Search" (sic).

Beberapa tahun yang lalu saya mengirim ini ke AWK. Versi ini mungkin lebih cocok untuk pengembang modern untuk port karena sintaksinya seperti C (daripada Fortran) dan beberapa komentar tambahan yang saya masukkan ketika mengembangkan dan mengujinya, karena saya perlu meningkatkan akurasinya. Itu muncul di bawah.

Bagi mereka yang tidak memiliki banyak pengalaman mengaitkan kode ilmiah / matematika / statistik, beberapa kata nasihat : satu kesalahan ketik dapat menyebabkan kesalahan serius yang mungkin tidak mudah terdeteksi. (Percayalah pada saya, saya sudah membuat banyak dari mereka.) Selalu, selalu membuat tes yang cermat dan lengkap. Karena fungsi integral / Gaussian terpisahkan / kesalahan normal tersedia dalam begitu banyak tabel dan begitu banyak perangkat lunak, mudah dan cepat untuk mentabulasikan sejumlah besar nilai fungsi porting Anda dan membandingkan secara sistematis (mis. Dengan komputer, bukan dengan mata) nilai untuk memperbaikinya. Anda dapat melihat tes seperti itu di awal kode saya: menghasilkan tabel nilai di -8.5: 8.5 (sebesar 0,1) yang dapat disalurkan (melalui STDOUT) ke program lain untuk pemeriksaan sistematis.

Pendekatan pengujian lain - bagi mereka yang memiliki latar belakang analisis numerik yang cukup untuk mengetahui bagaimana memperkirakan kesalahan yang diharapkan - akan secara numerik membedakan nilai-nilai dan membandingkannya dengan PDF (yang siap dihitung).

By the way: kode ini hanya untuk kasus dengan rata-rata $0$dan deviasi standar satuan ("sigma"). Tapi hanya itu yang dibutuhkan: untuk diintegrasikan$-\infty$ untuk $x$ ketika rata-rata adalah $\mu$ dan SD $\sigma$, hitung saja $z = (x-\mu)/\sigma$dan berlaku alnormuntuk itu.

Edit

Saya menguji port alnormto Mathematica, yang menghitung nilai-nilai dengan presisi sewenang-wenang. Untuk membandingkan hasilnya, berikut adalah plot log natural dari rasio nilai ekor atas$1 - \Phi(z)$ dengan $z\ge 1$. (Kesalahan relatif positif berarti alnormterlalu besar.)

Nilai selalu akurat untuk$4 \times 10^{-11}$ relatif terhadap probabilitas ekor yang semakin kecil . Anda dapat melihat di mana perhitungan beralih ke rumus asimptotik (di$z=16$) dan terbukti bahwa formula ini menjadi sangat akurat $z$meningkat. Plot berhenti di$z=\sqrt(2\times 708) \approx 37.6$ karena di sinilah eksponensial presisi ganda mulai kurang mengalir.

Misalnya, alnorm[-6.0]kembali$9.865\ 876\ 450\ 315E-10$ sedangkan nilai sebenarnya, sama dengan $\frac{1}{2}\text{erfc}(3\sqrt{2})$, kira-kira $9.865\ 876\ 450\ 377E-10$, pertama berbeda dalam digit desimal kedua belas.

NB Sebagai bagian dari pengeditan ini, saya mengubah UPPER_TAIL_IS_ZEROdari 15.ke 16.dalam kode: itu membuat hasilnya sedikit lebih akurat$Z$ antara $15$ dan $16$. (Akhir edit.)

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