¿Cómo calcular la probabilidad asociada con puntuaciones Z absurdamente grandes?

Question

Paquetes de Software para la red de detección motivo puede volver altísimas puntuaciones Z (el más alto que he visto es de 600.000+, pero las puntuaciones Z de más de 100 son bastante comunes). I plan para demostrar que estas puntuaciones Z son falsos.

Enorme puntuaciones Z corresponden a muy bajas probabilidades asociadas. Los valores de las probabilidades se dan, por ejemplo, en la distribución normal, en la página de la wikipedia (y probablemente todos los stats de libros de texto) para las puntuaciones Z de hasta 6. Así que...

Pregunta: ¿Cómo se puede calcular la función de error $1-\mathrm{erf}(n/\sqrt{2})$ n de hasta 1.000.000 de, por ejemplo?

Estoy particularmente después de una ya implementado de este paquete (si es posible). La mejor que he encontrado hasta ahora es WolframAlpha, que se las arregla para calcular que para n=150 (aquí).

Answer 1

Puede aproximar con mucho más simple de las funciones - ver esta Wikipedia sección para obtener más información. La aproximación básica es que el $\textrm{erf}(x) \approx \textrm{sgn}(x)\sqrt{1 - \exp(-x^2 \frac{4/\pi + ax^2}{1+ax^2}})$

El artículo tiene un enlace incorrecto para esa sección. El PDF que se hace referencia se puede encontrar en Sergei Winitzki de archivos - o en este enlace.

Answer 2

Un simple límite superior

Para valores muy grandes de la argumentación en el cálculo de la cola superior probabilidad de que una normal, excelente límites existen y que son, probablemente, como bueno como se va a obtener de utilizar otros métodos con doble precisión de punto flotante. Para $z > 0$, vamos $$\renewcommand{\Pr}{\mathbb{P}}\newcommand{\rd}{\mathrm{d}} S(z) := \Pr(Z > z) = \int_z^\infty \varphi(z) \rd z \>,$$ donde $\varphi(z) = (2\pi)^{-1/2} e^{-z^2/2}$ es el estándar normal pdf. He usado la notación $S(z)$, en deferencia a la norma la notación en el análisis de supervivencia. En la ingeniería de contextos, ellos llaman a esto la función de la $Q$-función y se denota por a $Q(z)$.

A continuación, una muy simple, elemental límite superior es $$\newcommand{\Ub}{\hat{S}_u} \newcommand{\Sl}{\hat{S}_\ell} S(z) \leq \frac{\varphi(z)}{z} =: \Ub(z) \> ,$$ donde la notación de la mano derecha indica que este es un -límite superior de la estimación. Esta respuesta da una prueba de la envolvente.

Hay varios buenos complementarios límites inferiores así. Uno de los más práctica y más fácil de derivar es el obligado $$S(z) \geq \frac{z}{z^2+1} \varphi(z) =: \Sl(z) \> .$$ Existen al menos tres métodos para la obtención de esta obligado. Un boceto de uno de estos métodos se puede encontrar en esta respuesta a un relacionados con la pregunta.

Una foto

Abajo es un dibujo de los dos límites (en gris), junto con el real la función $S(z)$.

¿Qué tan buena es?

De la trama, parece que los límites se vuelven bastante apretado, incluso para moderadamente grande $z$. Podríamos preguntarnos cómo son y apretado ¿qué tipo de cuantitativa declaración en ese sentido se puede hacer.

Una medida útil de la opresión es un absoluto error relativo $$\newcommand{\err}{\mathcal{E}} \err(z) = \left|\frac{\Ub(z) - S(z)}{S(z)}\right| \>.$$ Esto le da a usted el error de proporción de la estimación.

Ahora, tenga en cuenta que, dado que todas las funciones son no negativos, mediante la delimitación de las propiedades de $\Su(z)$$\Sl(z)$, obtenemos $$\err(z) = \frac{\Ub(z) - S(z)}{S(z)} \leq \frac{\Ub(z) - \Sl(z)}{\Sl(z)} = z^{-2} \> ,$$ y así esto proporciona una prueba de que para $z \geq 10$ el-límite superior es correcta dentro del 1%, $z \geq 28$ es correcto dentro de 0,1% y $z \geq 100$ es correcto dentro de 0.01%.

De hecho, la forma simple de los límites proporciona una buena manera de comprobar en otros "aproximaciones". Si, en el cálculo numérico de los más complicados aproximaciones, obtenemos un valor fuera de estos límites, podemos simplemente "correcto" que toma el valor de, por ejemplo, el límite superior se proporcionan aquí.

Hay muchos refinamientos de estos límites. El Laplace límites mencionados aquí proporcionan una agradable secuencia de los límites superior e inferior en $S(z)$ de la forma $R(z) \varphi(z)$ donde $R(z)$ es una función racional.

Por último, aquí es otro algo relacionado con la pregunta y la respuesta.

Answer 3

La pregunta se refiere a la función complementaria de error

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

for "large" values of $x$ ($=n/\sqrt{2}$ in the original question)--that is, between 100 and 700,000 or so. (In practice, any value greater than about 6 should be considered "large," as we will see.) Note that because this will be used to compute p-values, there is little value in obtaining more than three significant (decimal) digits.

To begin, consider the approximation suggested by @Iterator,

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

where

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

Although this is an excellent approximation to the error function itself, it's a terrible approximation to $\textrm{erfc}$. However, there is a way to systematically fix that up.

For the p-values associated with such large values of $x$, we're interested in the relative error $f(x)/\textrm{erfc}(x) - 1$: we hope its absolute value would be less than 0.001 for three significant digits of precision. Unfortunately this expression is difficult to study for large $x$ due to underflows in double-precision computation. Here is one attempt, which plots the relative error versus $x$ for $0 \le x \le 5.8$:

The calculation becomes unstable once $x$ exceeds 5.3 or so and cannot deliver one significant digit past 5.8. This is no surprise: $\exp(-5.8^2) \aprox 10^{-14.6}$ is pushing the limits of double-precision arithmetic. Because there is no evidence that the relative error is going to be acceptably small for larger $x$, we need to do better.

Performing the calculation in extended arithmetic (with Mathematica) improves our picture of what's going on:

The error increases rapidly with $x$ and shows no signs of leveling off. Past $x=10$ or so, this approximation doesn't even deliver one reliable digit of information!

However, the plot is starting to look linear. We might guess that the relative error is directly proportional to $x$. (This makes sense on theoretical grounds: $\textrm{erfc}$ is manifestly an odd function and $f$ is manifestly even, so their ratio ought to be an odd function. Thus we would expect the relative error, if it increases, to behave like an odd power of $x$.) This leads us to study the relative error divided by $x$. Equivalently, I choose to examine $x \cdot \textrm{erfc}(x)/f(x)$, because the hope is this should have a constant limiting value. Here is its graph:

Our guess appears to be borne out: this ratio does seem to be approaching a limit around 8 or so. When asked, Mathematica will supply it:

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

The value is $a_1 = \frac{2 }{\sqrt{\pi }}e^{\frac{3 (-4+\pi )^2}{8 (-3+\pi )}} \approx 7.94325$. This enables us to improve the estimate: we take

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

as the first refinement of the approximation. When $x$ is really large--greater than a few thousand--this approximation is just fine. Because it's still not going to be good enough for an interesting range of arguments between $5.3$ and $2000$ or so, let's iterate the procedure. This time, the inverse relative error--specifically, the expression $1 - \textrm{erfc}(x)/f_1(x)$--should behave like $1/x^2$ for large $x$ (by virtue of the previous parity considerations). Accordingly, we multiply by $x^2$ and find the next limit:

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

The value is

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

This process can proceed as long as we like. I took it out one more step, finding

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

with value approximately 1623.67. (The full expression involves a degree-eight rational function of $\pi$ and is too long to be useful here.)

Unwinding these operations yields our final approximation

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

The error is proportional to $x^{-6}$. Of import is the constant of proportionality, so we plot $x^6(1 - \textrm{erfc}(x) / f_3(x))$:

It rapidly approaches a limiting value around 2660.59. Using the approximation $f_3$, we obtain estimates of $\textrm{erfc}(x)$ whose relative accuracy is better than $2661/x^6$ for all $x \gt 0$. Once $x$ exceeds 20 or so, we have our three significant digits (or far more, as $x$ gets larger). As a check, here is a table comparing the correct values to the approximation for $x$ between $10$ and $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

In fact, this approximation delivers at least two significant figures of precision for $x=8$ on, which is just about where pedestrian calculations (such as Excel's NormSDist function) peter out.

Finally, one might worry about our ability to compute the initial approximation $f$. However, that's not hard: when $x$ is large enough to cause underflows in the exponential, the square root is well approximated by half the exponential,

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

Computing the logarithm of this (in base 10) is simple, and readily gives the desired result. For example, let $x=1000$. The common logarithm of this approximation is

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

Exponentiating yields

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

Applying the correction (in $f_3$) produces

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

Note that the correction reduces the original approximation by over 99% (and indeed, $a_1/x \aprox 1\text{%}$.) (This approximation differs from the correct value only in the last digit. Another well-known approximation, $\exp(-x^2)/(x\sqrt{\pi})$, equals $1.860038 \cdot 10^{-434298}$, errar en el sexto dígito significativo. Estoy seguro de que podemos mejorar, que también, si queremos, utilizando las mismas técnicas.)