R малые значения p

Я рассчитываю z-показатели, чтобы увидеть, далеко ли значение от среднего/медианы распределения. Первоначально я сделал это, используя среднее значение, а затем превратил их в двусторонние значения pvalue. Но теперь, используя медиану, я заметил, что в значениях p есть некоторые значения Na.

Я определил, что это происходит для значений, очень далеких от медианы. И похоже, что это связано с расчетом pnorm. «qnorm» основан на алгоритме Wichura AS 241, который обеспечивает точные результаты примерно до 16 цифр».

Кто-нибудь знает способ обойти это, так как мне нужны очень маленькие значения pvalue. Спасибо,

> z<- -12.5
> 2-2*pnorm(abs(z))
[1] 0
> z<- -10  
> 2-2*pnorm(abs(z))
[1] 0 
> z<- -8 
> 2-2*pnorm(abs(z))
[1] 1.332268e-15

person user2814482    schedule 20.10.2014    source источник


Ответы (3)


На самом деле вы вычисляете очень высокие p-значения:

options(digits=22)
z <- c(-12.5,-10,-8)
pnorm(abs(z))
# [1] 1.0000000000000000000000 1.0000000000000000000000 0.9999999999999993338662
2-2*pnorm(abs(z))
# [1] 0.000000000000000000000e+00 0.000000000000000000000e+00 1.332267629550187848508e-15

Я думаю, вам будет лучше использовать низкие p-значения (близкие к нулю), но я недостаточно хорош в математике, чтобы знать, связана ли ошибка с близкими к единице p-значениями в алгоритме AS241 или в хранилище с плавающей запятой. . Посмотрите, как красиво отображаются низкие значения:

pnorm(z)
# [1] 3.732564298877713761239e-36 7.619853024160526919908e-24 6.220960574271784860433e-16

Имейте в виду, что 1 - pnorm(x) эквивалентно pnorm(-x). Итак, 2-2*pnorm(abs(x)) эквивалентно 2*(1 - pnorm(abs(x)) эквивалентно 2*pnorm(-abs(x)), так что просто используйте:

2 * pnorm(-abs(z))
#  [1] 7.465128597755427522478e-36 1.523970604832105383982e-23 1.244192114854356972087e-15

который должен получить более точно то, что вы ищете.

person C8H10N4O2    schedule 20.10.2014

Одна мысль: вам придется использовать exp() с большей точностью, но вы могли бы использовать log(p), чтобы получить немного большую точность в хвостах, иначе вы фактически на 0 для нелогарифмических значений p с точки зрения диапазона, который может быть рассчитан:

> z<- -12.5
> pnorm(abs(z),log.p=T)
[1] -7.619853e-24

Преобразование обратно в значение p не работает, но вы можете сравнить его с log(p)...

> exp(pnorm(abs(z),log.p=T))
[1] 1
person Mike Burr    schedule 20.10.2014

pnorm — это функция, которая дает значение P на основе заданного x. Если вы не укажете больше аргументов, то распределение по умолчанию будет Нормальным со средним значением 0 и стандартным отклонением 1.
На основе симметрии pnorm(a) = 1-pnorm(-a).
Если в R добавить положительные числа, они будут округлены. Но если вы добавите минус, округление не будет выполнено. Таким образом, используя эту формулу и отрицательные числа, вы можете рассчитать необходимые значения.

> pnorm(0.25)
[1] 0.5987063
> 1-pnorm(-0.25)
[1] 0.5987063
> pnorm(20)
[1] 1
> pnorm(-20)
[1] 2.753624e-89
person K.I.    schedule 20.10.2014
comment
Мне кажется, что это не работает для сценария, описанного в вопросе 2-2*pnorm(abs(z)). - person SimonG; 20.10.2014