Bagaimana cara menghitung probabilitas yang terkait dengan skor-Z yang sangat besar?

Paket perangkat lunak untuk deteksi motif jaringan dapat mengembalikan skor Z yang sangat tinggi (nilai tertinggi yang pernah saya lihat adalah 600.000+, tetapi skor Z lebih dari 100 cukup umum). Saya berencana untuk menunjukkan bahwa skor-Z ini palsu.

Skor Z besar sesuai dengan probabilitas terkait yang sangat rendah. Nilai-nilai probabilitas terkait diberikan pada misalnya halaman wikipedia distribusi normal (dan mungkin setiap buku teks statistik) untuk skor-Z hingga 6. Jadi ...

Pertanyaan : Bagaimana cara menghitung fungsi kesalahan $1-\mathrm{erf}(n/\sqrt{2})$untuk n hingga 1.000.000, katakan?

Saya terutama setelah paket yang sudah diterapkan untuk ini (jika mungkin). Yang terbaik yang saya temukan sejauh ini adalah WolframAlpha, yang berhasil menghitungnya untuk n = 150 (di sini ).

Pertanyaannya menyangkut fungsi kesalahan komplementer

$$\textrm{erfc}(x) = \frac{2}{\sqrt{\pi}}\int_{x}^{\infty}\exp(-t^2) d t$$

untuk nilai "besar" ( = n / $x$ dalam pertanyaan asli) - yaitu, antara 100 dan 700.000 atau lebih. (Dalam praktiknya, nilai apa pun yang lebih besar dari sekitar 6 harus dianggap "besar," seperti yang akan kita lihat.) Perhatikan bahwa karena ini akan digunakan untuk menghitung nilai-p, ada sedikit nilai dalam memperoleh lebih dari tiga digit signifikan (desimal) .$=n/\sqrt{2}$

Untuk memulai, pertimbangkan perkiraan yang disarankan oleh @Iterator,

$$f(x) = 1 - \sqrt{1 - \exp \left(-x^2 \left(\frac{4 + ax^2}{\pi + ax^2}\right)\right)},$$

dimana

$$a = \frac{8(\pi-3)}{3(4-\pi)} \approx 0.439862.$$

Meskipun ini merupakan perkiraan yang sangat baik untuk fungsi kesalahan itu sendiri, ini merupakan perkiraan yang mengerikan untuk . Namun, ada cara untuk memperbaikinya secara sistematis.$\textrm{erfc}$

Untuk nilai-p yang terkait dengan nilai sedemikian besar , kami tertarik pada kesalahan relatif f ( x ) / erfc ( x ) - 1 : kami berharap nilai absolutnya akan kurang dari 0,001 untuk tiga digit presisi yang signifikan. Sayangnya ungkapan ini sulit dipelajari untuk x besar$x$ $f(x)/\textrm{erfc}(x) - 1$$x$ karena arus bawah dalam perhitungan presisi ganda. Berikut adalah satu upaya, yang menggambarkan kesalahan relatif versus untuk 0 ≤ x ≤ 5.8 :$x$$0 \le x \le 5.8$

Perhitungan menjadi tidak stabil setelah melebihi 5,3 atau lebih dan tidak dapat memberikan satu digit signifikan melewati 5,8. Ini tidak mengherankan: exp ( - $x$ mendorong batas-batas aritmatika presisi ganda. Karena tidak ada bukti bahwa kesalahan relatif akan menjadi kecil untuk x yang lebih besar, kita perlu melakukan yang lebih baik.$\exp(-5.8^2) \approx 10^{-14.6}$$x$

Melakukan perhitungan dalam aritmatika diperpanjang (dengan Mathematica ) meningkatkan gambaran kita tentang apa yang terjadi:

Kesalahan meningkat dengan cepat dengan dan tidak menunjukkan tanda-tanda leveling off. Melewati x = 10 atau lebih, perkiraan ini bahkan tidak memberikan satu digit informasi yang dapat diandalkan!$x$$x=10$

Namun, alurnya mulai terlihat linier. Kami mungkin menduga bahwa kesalahan relatif berbanding lurus dengan . (Ini masuk akal dengan alasan teoritis: erfc secara nyata merupakan fungsi yang ganjil dan f adalah bahkan secara nyata, sehingga rasio mereka seharusnya merupakan fungsi yang ganjil. Dengan demikian kita akan mengharapkan kesalahan relatif, jika meningkat, untuk berperilaku seperti kekuatan aneh x .) Ini menuntun kita untuk mempelajari kesalahan relatif dibagi x . Setara, saya memilih untuk memeriksa x ⋅ erfc ( x ) / f ( x $x$$\textrm{erfc}$$f$$x$ $x$$x \cdot \textrm{erfc}(x)/f(x)$, karena harapannya adalah ini harus memiliki nilai pembatas yang konstan. Ini grafiknya:

Dugaan kami tampaknya terbukti: rasio ini tampaknya mendekati batas sekitar 8 atau lebih. Ketika ditanya, Mathematica akan menyediakannya:

a1 = Limit[x (Erfc[x]/f[x]), x -> \[Infinity]]

Nilainya adalah . Ini memungkinkan kami untuk meningkatkan taksiran:kami ambil$a_1 = \frac{2 }{\sqrt{\pi }}e^{\frac{3 (-4+\pi )^2}{8 (-3+\pi )}} \approx 7.94325$

$$f_1(x) = f(x) \frac{a_1}{x}$$

sebagai penyempurnaan pertama aproksimasi. Ketika benar-benar besar - lebih besar dari beberapa ribu - perkiraan ini baik-baik saja. Karena itu masih tidak akan cukup baik untuk berbagai argumen yang menarik antara 5,3 dan 2000 atau lebih, mari kita beralih prosedur. Kali ini, kesalahan relatif terbalik - khususnya, ekspresi 1 - erfc ( x ) / f 1 ( x ) - harus berperilaku seperti 1 / x 2 untuk x besar (berdasarkan pertimbangan paritas sebelumnya). Dengan demikian, kita kalikan dengan x $x$$5.3$$2000$$1 - \textrm{erfc}(x)/f_1(x)$$1/x^2$$x$$x^2$ dan temukan batas selanjutnya:

a2 = Limit[x^2 (a1 - x (Erfc[x]/f[x])), x -> \[Infinity]] 

Nilainya adalah

$$a_2 = \frac{1}{32 \sqrt{\pi }} e^{\frac{3 (-4+\pi )^2}{8 (-3+\pi )}} \left(32-\frac{9 (-4+\pi )^3 \pi }{(-3+\pi )^2}\right) \approx 114.687.$$

Proses ini dapat berlangsung selama yang kita mau. Saya mengambilnya satu langkah lagi, menemukan

a3 = Limit[x^2 (a2 - x^2 (a1 - x (Erfc[x]/f[x]))), x -> \[Infinity]] 

dengan nilai sekitar 1623.67. (Ekspresi penuh melibatkan fungsi rasional tingkat delapan dari dan terlalu lama untuk berguna di sini.)$\pi$

Mengurai operasi ini menghasilkan perkiraan akhir kami

$$f_3(x) = f(x)\left(a_1 - a_2/x^2 + a_3/x^4\right)/x.$$

Kesalahan sebanding dengan . Impor adalah konstanta proporsionalitas, jadi kami memplot x 6 ( 1 - erfc ( x ) / f 3 ( x ) ) :$x^{-6}$$x^6(1 - \textrm{erfc}(x) / f_3(x))$

Dengan cepat mendekati nilai pembatas sekitar 2660.59. Dengan menggunakan perkiraan , kami memperoleh estimasi erfc ( x ) yang akurasi relatifnya lebih baik dari 2661 / x 6 untuk semua x > 0 . Setelah x melebihi 20 atau lebih, kita memiliki tiga digit signifikan (atau lebih, karena x semakin besar). Sebagai tanda centang, berikut adalah tabel yang membandingkan nilai yang benar dengan perkiraan untuk x antara 10 dan $f_3$$\textrm{erfc}(x)$$2661/x^6$$x \gt 0$$x$$x$$x$$10$$20$ :

 x  Erfc    Approximation      
10  2.088*10^-45    2.094*10^-45
11  1.441*10^-54    1.443*10^-54
12  1.356*10^-64    1.357*10^-64
13  1.740*10^-75    1.741*10^-75
14  3.037*10^-87    3.038*10^-87
15  7.213*10^-100   7.215*10^-100
16  2.328*10^-113   2.329*10^-113
17  1.021*10^-127   1.021*10^-127
18  6.082*10^-143   6.083*10^-143
19  4.918*10^-159   4.918*10^-159
20  5.396*10^-176   5.396*10^-176

Bahkan, perkiraan ini memberikan setidaknya dua angka presisi yang signifikan untuk pada, yang hanya tentang di mana perhitungan pejalan kaki (sepertifungsiExcel) peter out.$x=8$NormSDist

Akhirnya, orang mungkin khawatir tentang kemampuan kita untuk menghitung perkiraan awal . Namun, itu tidak sulit: ketika x cukup besar untuk menyebabkan arus bawah dalam eksponensial, akar kuadrat diperkirakan hampir separuh eksponensial,$f$$x$

$$f(x) \approx \frac{1}{2} \exp \left(-x^2 \left(\frac{4 + ax^2}{\pi + ax^2}\right)\right).$$

$x=1000$ . Logaritma umum dari pendekatan ini adalah

$$\log_{10}(f(x)) \approx \left(-1000^2 \left(\frac{4 + a \cdot 1000^2}{\pi + a \cdot 1000^2}\right) - \log(2)\right) / \log(10) \sim -434295.63047.$$

Hasil yang eksponensial

$$f(1000) \approx 2.34169 \cdot 10^{-434296}.$$

$f_3$

$$\textrm{erfc}(1000) \approx 1.86003\ 70486\ 32328 \cdot 10^{-434298}.$$

$a_1/x \approx 1\text{%}$$\exp(-x^2)/(x\sqrt{\pi})$$1.860038 \cdot 10^{-434298}$

Batas atas yang sederhana

$z > 0$

$$\renewcommand{\Pr}{\mathbb{P}}\newcommand{\rd}{\mathrm{d}} S(z) := \Pr(Z > z) = \int_z^\infty \varphi(z) \rd z \>,$$ $\varphi(z) = (2\pi)^{-1/2} e^{-z^2/2}$$S(z)$$Q$$Q(z)$

$$\newcommand{\Su}{\hat{S}_u} \newcommand{\Sl}{\hat{S}_\ell} S(z) \leq \frac{\varphi(z)}{z} =: \Su(z) \> ,$$

$$S(z) \geq \frac{z}{z^2+1} \varphi(z) =: \Sl(z) \> .$$

Sebuah gambar

$S(z)$

Seberapa baik itu?

$z$

$$\newcommand{\err}{\mathcal{E}} \err(z) = \left|\frac{\Su(z) - S(z)}{S(z)}\right| \>.$$ Ini memberi Anda kesalahan proporsional estimasi.

$\Su(z)$$\Sl(z)$

$$\err(z) = \frac{\Su(z) - S(z)}{S(z)} \leq \frac{\Su(z) - \Sl(z)}{\Sl(z)} = z^{-2} \> ,$$ $z \geq 10$$z \geq 28$$z \geq 100$ itu benar dalam 0,01%.

Bahkan, bentuk batas yang sederhana memberikan pemeriksaan yang baik pada "perkiraan" lainnya. Jika, dalam perhitungan numerik perkiraan yang lebih rumit, kita mendapatkan nilai di luar batas ini, kita bisa "mengoreksi" itu untuk mengambil nilai, misalnya, batas atas yang disediakan di sini.

$S(z)$$R(z) \varphi(z)$$R(z)$ adalah fungsi rasional.

Akhirnya, inilah pertanyaan dan jawaban yang agak terkait.

Anda dapat memperkirakannya dengan fungsi yang lebih sederhana - lihat bagian Wikipedia ini untuk informasi lebih lanjut. Perkiraan dasarnya adalah itu$\textrm{erf}(x) \approx \textrm{sgn}(x)\sqrt{1 - \exp(-x^2 \frac{4/\pi + ax^2}{1+ax^2}})$

Artikel memiliki tautan yang salah untuk bagian itu. Referensi PDF dapat ditemukan dalam file Sergei Winitzki - atau di tautan ini .