Как рассчитать вероятность, связанную с нелепо большими Z-показателями?

14

Пакеты программ для обнаружения сетевых мотивов могут возвращать чрезвычайно высокие Z-оценки (самый высокий показатель, который я видел, составляет 600 000+, но Z-оценки более 100 встречаются довольно часто). Я планирую показать, что эти Z-оценки являются поддельными.

Огромные Z-оценки соответствуют чрезвычайно низким вероятностям. Значения связанных вероятностей приведены, например, на странице википедии нормального распределения (и, вероятно, в каждом учебнике статистики) для Z-показателей до 6. Итак ...

Вопрос : Как вычислить функцию ошибки 1erf(n/2)для n до 1000000, говорите?

Я особенно после уже реализованного пакета для этого (если это возможно). Лучшее, что я нашел до сих пор, это WolframAlpha, которому удается вычислить его для n = 150 ( здесь ).

Дуглас С. Стоунс
источник
6
Может быть, это не тот вопрос, который нужно задавать. Эти z-показатели являются поддельными, поскольку они предполагают, что нормальное распределение является гораздо лучшим приближением или моделью, чем оно есть на самом деле. Это немного похоже на предположение, что ньютоновская механика хороша до 600 000 десятичных знаков. Если вы действительно заинтересованы исключительно в вычислении erf для экстремальных значений n , то этот вопрос относится к math.SE, а не здесь.
whuber
6
Для «нелепо» больших значений вы не добьетесь большего успеха, чем использование верхней границы Pr(Z>z)(z2π)1ez2/2 для двойной точностиплавающей точкой. Это приближение и другие обсуждаются в другом месте на stats.SE.
кардинал
Спасибо, кардинал, эта оценка кажется довольно точной. Почему вы не делаете это ответ?
Дуглас С. Стоунс
@Douglas: Если вы все еще заинтересованы, я могу что-то собрать на следующий день или около того и опубликовать это как более полный ответ.
кардинал
1
Ну ... я думаю, что стоит добавить это как ответ. Может быть, границы - это общеизвестно в prob + stats, но я этого не знал. Кроме того, Q и A здесь не только для OP.
Дуглас С. Стоунс

Ответы:

19

Вопрос касается дополнительной функции ошибок

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

для «больших» значений ( = n / x в исходном вопросе) - то есть между 100 и 700 000 или около того. (На практике любое значение, большее, чем примерно 6, следует рассматривать как «большое», как мы увидим.) Обратите внимание, что поскольку оно будет использоваться для вычисления p-значений, при получении более трех значащих (десятичных) цифр мало значения. ,=n/2

Для начала рассмотрим приближение, предложенное @Iterator,

f(x)=11exp(x2(4+ax2π+ax2)),

где

a=8(π3)3(4π)0.439862.

Хотя это отличное приближение к самой функции ошибки, это ужасное приближение к . Однако есть способ систематически исправить это.erfc

Для значений p, связанных с такими большими значениями , нас интересует относительная ошибка f ( x ) / erfc ( x ) - 1 : мы надеемся, что ее абсолютное значение будет менее 0,001 для трех значащих цифр точности. К сожалению, это выражение трудно изучить для больших x из-за недостаточного вычисления при двойной точности. Вот одна попытка, которая показывает относительную ошибку по сравнению с x для 0 x 5,8 :x f(x)/erfc(x)1xx0x5.8

Участок 1

Расчет становится нестабильным, когда превышает 5,3 или около того и не может поставить одну значащую цифру после 5,8. Это неудивительно: exp ( - 5.8 2 ) 10 - 14.6 расширяет границы арифметики двойной точности. Поскольку нет никаких доказательств того, что относительная погрешность будет приемлемо мала для больших x , мы должны сделать лучше.xexp(5.82)1014.6x

Выполнение вычислений в расширенной арифметике (с Mathematica ) улучшает нашу картину происходящего:

Участок 2

Ошибка быстро увеличивается с и не показывает никаких признаков выравнивания. В прошлом х = 10 или около того, это приближение даже не доставляет одну достоверную цифру информации!xx=10

Тем не менее, сюжет начинает выглядеть линейно. Можно предположить, что относительная ошибка прямо пропорциональна . (Это имеет смысл с теоретической точки зрения : erfc является явно нечетной функцией, а f явно четной, поэтому их отношение должно быть нечетной функцией. Таким образом, можно ожидать, что относительная ошибка, если она возрастет, будет вести себя как нечетная степень x .) Это приводит нас к изучению относительной ошибки, деленной на х . Эквивалентно, я выбираю изучить x erfc ( x ) / f ( x )xerfcfx xxerfc(x)/f(x)Потому что надежда на это должна иметь постоянное предельное значение. Вот его график:

Участок 3

Наше предположение, похоже, подтверждается: это соотношение, похоже, приближается к пределу около 8 или около того. Когда спрошено, Mathematica предоставит это:

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

Значение 1 = 2. Это позволяет нам улучшить оценку:мы беремa1=2πe3(4+π)28(3+π)7.94325

f1(x)=f(x)a1x

в качестве первого уточнения приближения. Когда действительно велико - больше нескольких тысяч - это приближение просто отлично. Поскольку он все еще не будет достаточно хорош для интересного диапазона аргументов между 5.3 и 2000 или около того, давайте повторим процедуру. На этот раз обратная относительная ошибка - в частности, выражение 1 - erfc ( x ) / f 1 ( x ) - должна вести себя как 1 / x 2 для больших x (в силу предыдущих соображений четности). Соответственно умножаем на х 2x5.320001erfc(x)/f1(x)1/x2xx2 и найдите следующий предел:

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

Значение

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

Этот процесс может продолжаться сколько угодно. Я сделал еще один шаг, найдя

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

со значением приблизительно 1623,67. (Полное выражение включает рациональную функцию степени восьмого и слишком длинное, чтобы быть здесь полезным.)π

Разматывание этих операций дает наше окончательное приближение

f3(x)=f(x)(a1a2/x2+a3/Икс4)/Икс,

The error is proportional to Икс-6. Of import is the constant of proportionality, so we plot Икс6(1-ERFC(Икс)/е3(Икс)):

Сюжет 4

It rapidly approaches a limiting value around 2660.59. Using the approximation е3, we obtain estimates of ERFC(Икс) whose relative accuracy is better than 2661/Икс6 for all Икс>0. Once Икс exceeds 20 or so, we have our three significant digits (or far more, as Икс gets larger). As a check, here is a table comparing the correct values to the approximation for Икс 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

Иксзнак равно8NormSDist

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)12exp(x2(4+ax2π+ax2)).

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

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
источник
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).

Верхний хвост нормалей и границ

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
источник
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)

Статья содержит неправильную ссылку на этот раздел. Ссылочный PDF-файл можно найти в файлах Сергея Виницкого или по этой ссылке .

Итератор
источник
1
Некоторое усиление этого было бы желательно по двум причинам. Во-первых, лучше, когда ответы могут стоять в одиночестве. Во-вторых, эта статья неоднозначно пишет о качестве аппроксимации «в бесконечности»: насколько точной является «очень точная»? (У вас явно есть это хорошее понимание, но от всех заинтересованных читателей этого ожидать не стоит.) Заявленное значение «.00035» здесь бесполезно.
Whuber
Благодарю. Я не заметил, что была поддержка на основе Javascript для использования TeX, что имело значение при написании этого.
Итератор
1
Кстати, ссылка на Википедию в этом приближении не работает. Однако Mathematica считает, что относительная ошибка (1 - приблизительно (x) / erf (x)) ведет себя как обратная величина2exp(Икс2+3(π-4)2/(8(π-3))),
whuber
@whuber, можете ли вы опубликовать код Mathematica для этого? :) Я не видел Mathematica более 15 лет, и никогда для такого рода целей.
Итератор
Я разместил это в отдельном ответе.
whuber