Como calcular a probabilidade associada a escores Z absurdamente grandes?

14

Pacotes de software para detecção de motivo de rede podem retornar escores Z enormemente altos (o mais alto que já vi são mais de 600.000, mas escores Z de mais de 100 são bastante comuns). Eu pretendo mostrar que esses Z-scores são falsos.

Escores Z enormes correspondem a probabilidades associadas extremamente baixas. Os valores das probabilidades associadas são dados, por exemplo, na página da Wikipedia de distribuição normal (e provavelmente em todos os livros de estatística) para pontuações Z de até 6. Então ...

Pergunta : Como calcular a função de erro 1erf(n/2)para n até 1.000.000, digamos?

Estou particularmente após um pacote já implementado para isso (se possível). O melhor que encontrei até agora é o WolframAlpha, que consegue calculá-lo para n = 150 ( aqui ).

Douglas S. Stones
fonte
6
Talvez essa não seja a pergunta certa a ser feita. Esses escores z são falsos porque assumem que a distribuição normal é uma aproximação ou modelo muito melhor do que realmente é. É como assumir que a mecânica newtoniana é boa para 600.000 casas decimais. Se você está realmente interessado apenas em calcular erf para valores extremos de , essa pergunta pertence a math.SE, não aqui. n
whuber
6
Para valores "absurdamente" grandes, você não fará melhor do que usar o limite superior para o ponto flutuante de precisão dupla. Essa aproximação e outras são discutidas em outras partes no stats.SE. Pr(Z>z)(z2π)1ez2/2
cardeal
Obrigado cardeal, esse limite parece ser bastante preciso. Por que você não faz disso uma resposta?
Douglas S. Stones
@ Douglas: Se você ainda estiver interessado, posso montar algo no dia seguinte ou publicá-lo como uma resposta mais completa.
cardeal
1
Bem ... acho que valeria a pena adicioná-lo como resposta. Talvez o limite seja o conhecimento comum em prob + stats, mas eu não sabia. Além disso, os Q e A aqui não são apenas para o OP.
Douglas S. Stones

Respostas:

19

A questão diz respeito à função de erro complementar

erfc(x)=2πxexp(-t2)dt

para valores "grandes" de (x na pergunta original) - ou seja, entre 100 e 700.000 ou mais. (Na prática, qualquer valor maior que cerca de 6 deve ser considerado "grande", como veremos.) Observe que, como isso será usado para calcular valores de p, há pouco valor na obtenção de mais de três dígitos significativos (decimais) .=n/2

Para começar, considere a aproximação sugerida pelo @Iterator,

f(x)=1-1-exp(-x2(4+umax2π+umax2)),

Onde

uma=8(π-3)3(4-π)0,439862.

Embora esta seja uma excelente aproximação à própria função de erro, é uma terrível aproximação à . No entanto, existe uma maneira de corrigir isso sistematicamente.erfc

Para os valores p associados a valores tão grandes de , estamos interessados ​​no erro relativo f ( x ) / erfc ( x ) - 1 : esperamos que seu valor absoluto seja menor que 0,001 para três dígitos significativos de precisão. Infelizmente, essa expressão é difícil de estudar para x grande devido a fluxos insuficientes no cálculo de precisão dupla. Aqui está uma tentativa, que plota o erro relativo versus x para 0 x x f(x)/erfc(x)-1xx :0 0x5,8

Gráfico 1

O cálculo se torna instável quando excede 5,3 ou mais e não pode fornecer um dígito significativo além de 5,8. Isso não é surpresa: exp ( - 5,8 2 ) 10 - 14,6 está ultrapassando os limites da aritmética de precisão dupla. Como não há evidências de que o erro relativo seja aceitável pequeno para um x maior , precisamos fazer melhor.xexp(-5,82)10-14,6x

A realização do cálculo em aritmética estendida (com o Mathematica ) melhora nossa imagem do que está acontecendo:

Gráfico 2

O erro aumenta rapidamente com e não mostra sinais de nivelamento. Passado xx , aproximadamente, essa aproximação nem fornece um dígito confiável de informações!x=10

No entanto, o enredo está começando a parecer linear. Podemos supor que o erro relativo seja diretamente proporcional a . (Isso faz sentido em bases teóricas: erfc é manifestamente uma função ímpar ef é manifestamente uniforme, portanto a razão deve ser uma função ímpar. Assim, esperaríamos que o erro relativo, se aumentasse, se comportasse como uma potência ímpar de x .) Isso nos leva a estudar o erro relativo dividido por x . Equivalentemente, eu escolho para examinar x ERFC ( x ) / f ( , porque a esperança é este deve ter um valor limite constante Aqui é o seu gráfico.:xerfcfx xxerfc(x)/f(x)

Gráfico 3

Nosso palpite parece estar confirmado: essa proporção parece estar se aproximando de um limite em torno de 8 ou mais. Quando solicitado, o Mathematica fornecerá:

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

O valor é . Isso nos permite melhorar a estimativa:tomamosa1=2πe3(4+π)28(3+π)7.94325

f1(x)=f(x)a1x

como o primeiro refinamento da aproximação. Quando é realmente grande - superior a alguns milhares - essa aproximação é ótima. Como ainda não será bom o suficiente para uma gama interessante de argumentos entre 5.3 e 2000 , vamos repetir o procedimento. Desta vez, o erro relativo inverso - especificamente, a expressão 1 - erfc ( x ) / f 1 ( x ) - deve se comportar como 1 / x 2 para x grande (em virtude das considerações de paridade anteriores). Assim, multiplicamos por x 2x5.320001erfc(x)/f1(x)1/x2xx2 e encontre o próximo limite:

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

O valor é

a2=132πe3(4+π)28(3+π)(329(4+π)3π(3+π)2)114.687.

Esse processo pode prosseguir o quanto quisermos. Eu dei mais um passo, encontrando

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

com valor aproximadamente 1623,67. (A expressão completa envolve uma função racional de grau oito de e é muito longa para ser útil aqui.)π

Desenrolar essas operações produz nossa aproximação final

f3(x)=f(x)(a1a2/x2+a3/x4)/x.

O erro é proporcional a . De importação é a constante de proporcionalidade, portanto, plotamos x 6 ( 1 - erfc ( x ) / f 3 ( x ) ) :x6x6(1erfc(x)/f3(x))

Gráfico 4

Ele se aproxima rapidamente de um valor limite em torno de 2660,59. Usando a aproximação , obtemos estimativas de erfc ( x ) cuja precisão relativa é melhor que 2661 / x 6 para todos x > 0 . Uma vez que x é superior a 20 ou assim, nós temos os nossos três dígitos significativos (ou muito mais, como x fica maior). Como verificação, segue uma tabela comparando os valores corretos com a aproximação de x entre 10 e 20 :f3erfc(x)2661/x6x>0xxx1020

 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

De fato, essa aproximação fornece pelo menos duas figuras significativas de precisão para x=8 , que é exatamente onde os cálculos de pedestres (como a NormSDistfunção do Excel ) desaparecem.

Finalmente, pode-se preocupar com nossa capacidade de calcular a aproximação inicial . No entanto, isso não é difícil: quando x é grande o suficiente para causar subfluxos no exponencial, a raiz quadrada é bem aproximada pela metade do exponencial,fx

f(x)12exp(x2(4+ax2π+ax2)).

A computação do logaritmo disso (na base 10) é simples e fornece rapidamente o resultado desejado. Por exemplo, deixe . O logaritmo comum dessa aproximação éx=1000

log10(f(x))(10002(4+a10002π+a10002)log(2))/log(10)434295.63047.

Exponentiating yields

f(1000)2.3416910434296.

Applying the correction (in f3) produces

erfc(1000)1.86003 70486 3232810434298.

Note that the correction reduces the original approximation by over 99% (and indeed, a1/x1%.) (This approximation differs from the correct value only in the last digit. Another well-known approximation, exp(x2)/(xπ), equals 1.86003810434298, erring in the sixth significant digit. I'm sure we could improve that one, too, if we wanted, using the same techniques.)

whuber
fonte
1
+1 This is a great answer, somehow I have never come across this thread before.
amoeba says Reinstate Monica
15

A simple upper bound

For very large values of the argument in the calculation of upper tail probability of a normal, excellent bounds exist that are probably as good as one will get using any other methods with double-precision floating point. For z>0, let

S(z):=P(Z>z)=zφ(z)dz,
where φ(z)=(2π)1/2ez2/2 is the standard normal pdf. I've used the notation S(z) in deference to the standard notation in survival analysis. In engineering contexts, they call this function the Q-function and denote it by Q(z).

Then, a very simple, elementary upper bound is

S(z)φ(z)z=:S^u(z),
where the notation on the right-hand side indicates this is an upper-bound estimate. This answer gives a proof of the bound.

There are several nice complementary lower bounds as well. One of the handiest and easiest to derive is the bound

S(z)zz2+1φ(z)=:S^(z).
There are at least three separate methods for deriving this bound. A rough sketch of one such method can be found in this answer to a related question.

A picture

Below is a plot of the two bounds (in grey) along with the actual function S(z).

Upper-tail of normal and bounds

How good is it?

From the plot, it seems that the bounds become quite tight even for moderately large z. We might ask ourselves how tight they are and what sort of quantitative statement in that regard can be made.

One useful measure of tightness is the absolute relative error

E(z)=|S^u(z)S(z)S(z)|.
This gives you the proportional error of the estimate.

Now, note that, since all of the involved functions are nonnegative, by using the bounding properties of S^u(z) and S^(z), we get

E(z)=S^u(z)S(z)S(z)S^u(z)S^(z)S^(z)=z2,
and so this provides a proof that for z10 the upper-bound is correct to within 1%, for z28 it is correct to within 0.1% and for z100 it is correct to within 0.01%.

In fact, the simple form of the bounds provides a good check on other "approximations". If, in the numerical calculation of more complicated approximations, we get a value outside these bounds, we can simply "correct" it to take the value of, e.g., the upper bound provided here.

There are many refinements of these bounds. The Laplace bounds mentioned here provide a nice sequence of upper and lower bounds on S(z) of the form R(z)φ(z) where R(z) is a rational function.

Finally, here is another somewhat-related question and answer.

cardinal
fonte
1
Apologies for all the "self-citations". Once, several years ago, I took an intense, two-week-long interest in related questions and tried to learn as much as I could about this topic.
cardinal
+1 Agree with whuber. Very nice, and I appreciate the links to other answers.
Iterator
5

You can approximate it with much simpler functions - see this Wikipedia section for more information. The basic approximation is that erf(x)sgn(x)1exp(x24/π+ax21+ax2)

The article has an incorrect link for that section. The PDF referenced can be found in Sergei Winitzki's files - or at this link.

Iterator
fonte
1
Some amplification of this would be welcome, for two reasons. First, it's best when answers can stand alone. Second, that article writes ambiguously about the quality of the approximation "in a neighborhood of infinity": just how accurate is "very accurate"? (You implicitly have a good sense of this, but it's a lot to expect of all interested readers.) The stated value of ".00035" is useless here.
whuber
Thanks. I didn't notice that there was Javascript-based support for using TeX, which made the difference in writing that out.
Iterator
1
Incidentally, the Wikipedia reference to that approximation is broken. Mathematica finds, though, that the relative error (1 - approx(x)/erf(x)) behaves like the reciprocal of 2exp(x2+3(π4)2/(8(π3))).
whuber
@whuber, você pode postar o código do Mathematica para isso? :) Não vejo o Mathematica há mais de 15 anos e nunca para esse tipo de objetivo.
Iterator
I posted it in a separate reply.
whuber