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 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 ).
fonte
Respostas:
A questão diz respeito à função de erro complementar
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,
Onde
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)−1 x x :0 ≤ x ≤ 5,8
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.x exp( - 5,82) ≈ 10- 14,6 x
A realização do cálculo em aritmética estendida (com o Mathematica ) melhora nossa imagem do que está acontecendo:
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.:x erfc f x x x⋅erfc(x)/f(x)
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á:
O valor é . Isso nos permite melhorar a estimativa:tomamosa1=2π√e3(−4+π)28(−3+π)≈7.94325
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 2x 5.3 2000 1−erfc(x)/f1(x) 1/x2 x x2 e encontre o próximo limite:
O valor é
Esse processo pode prosseguir o quanto quisermos. Eu dei mais um passo, encontrando
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
O erro é proporcional a . De importação é a constante de proporcionalidade, portanto, plotamos x 6 ( 1 - erfc ( x ) / f 3 ( x ) ) :x−6 x6(1−erfc(x)/f3(x))
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 :f3 erfc(x) 2661/x6 x>0 x x x 10 20
De fato, essa aproximação fornece pelo menos duas figuras significativas de precisão parax=8 , que é exatamente onde os cálculos de pedestres (como a
NormSDist
funçã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,f x
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
Exponentiating yields
Applying the correction (inf3 ) produces
Note that the correction reduces the original approximation by over 99% (and indeed,a1/x≈1% .) (This approximation differs from the correct value only in the last digit. Another well-known approximation, exp(−x2)/(xπ−−√) , equals 1.860038⋅10−434298 , erring in the sixth significant digit. I'm sure we could improve that one, too, if we wanted, using the same techniques.)
fonte
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. Forz>0 , let
Then, a very simple, elementary upper bound is
There are several nice complementary lower bounds as well. One of the handiest and easiest to derive is the bound
A picture
Below is a plot of the two bounds (in grey) along with the actual functionS(z) .
How good is it?
From the plot, it seems that the bounds become quite tight even for moderately largez . 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
Now, note that, since all of the involved functions are nonnegative, by using the bounding properties ofS^u(z) and S^ℓ(z) , we get
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 onS(z) of the form R(z)φ(z) where R(z) is a rational function.
Finally, here is another somewhat-related question and answer.
fonte
You can approximate it with much simpler functions - see this Wikipedia section for more information. The basic approximation is thaterf(x)≈sgn(x)1−exp(−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.
fonte