Использование круглой функции в R для очень маленьких значений возвращает ноль

Иногда мне приходится иметь дело с очень низкими значениями p и представлять их в табличном формате. Значения, возвращаемые R, могут иметь длинные значащие цифры (т. е. цифры после запятой). Теперь, поскольку значение p в любом случае такое низкое, я обычно сокращаю их, прежде чем записывать в файлы .xls или .tsv. (Просто чтобы таблицы выглядели красиво !!)

Я использую R version 3.0.0 (2013-04-03)

Немного предыстории и примеров:

9.881313e-208 в R будет 9.88e-208 в моей таблице

Для этого я могу использовать функцию round в R.

round(9.881313e-208, 210)
[1] 9.88e-208

Однако значение мощности e может различаться в каждом случае, и, поскольку таких случаев много, я использую следующую формулу:

x = 9.881313e-208
round(x,abs(floor(log10(x)-2)))   ## I came to this following trial and error
[1] 9.88e-208

Я проверил эту формулу эмпирически, и она работает в разных случаях, например:

a <- c(1.345678e-150,8.543678e-250,5.555555e-303, 0.01123, 4.523456e-290)
round(a,abs(floor(log10(a)-2)))
[1] 1.35e-150 8.54e-250 5.56e-303  1.12e-02 4.52e-290

Теперь проблема начинается, когда мощность e превышает число 306 (даже 307 нормально, но после 308 начинает становиться странно)

## Example 1:
b <- c(1.345678e-306,1.345678e-307,1.345678e-308, 1.345678e-309, 1.345678e-310)
round(b,abs(floor(log10(b)-2)))
[1] 1.35e-306 1.30e-307 1.00e-308  0.00e+00  0.00e+00

## Example 2:
b <- c(3.345678e-306,3.345678e-307,3.345678e-308, 3.345678e-309, 3.345678e-310)
round(b,abs(floor(log10(b)-2)))

[1] 3.35e-306 3.30e-307 3.00e-308  0.00e+00  0.00e+00

## Example 3:
b <- c(7.345678e-306,7.345678e-307,7.345678e-308, 7.345678e-309, 7.345678e-310)
round(b,abs(floor(log10(b)-2)))

[1] 7.35e-306 7.30e-307 7.00e-308 1.00e-308  0.00e+00

Кроме того, я проверил это напрямую:

round(7.356789e-306,308)
[1] 7.36e-306

round(7.356789e-307,309)
[1] 7.4e-307

round(7.356789e-308,310)
[1] 7e-308

round(7.356789e-309,311)
[1] 1e-308

round(7.356789e-310,312)
[1] 0

round(7.356789e-311,313)
[1] 0

Я пропустил что-то тривиальное здесь или функция round достигает предела разрешения за пределами e-308. Я знаю, что эти значения чрезвычайно малы и почти равны нулю, однако я бы все же предпочел иметь точное значение. Я видел некоторые ответы в SO с использованием Python для этой проблемы (см. Как правильно округлить очень маленькое отрицательное число?), но есть ли какие-нибудь предложения относительно того, как я могу преодолеть это в R?

Помощь очень ценится

Ваше здоровье

Эшвин


person Ashwin    schedule 21.05.2014    source источник
comment
Вы пробовали format() вместо round()?   -  person Andrie    schedule 21.05.2014
comment
Я согласен с Эндри: round следует использовать для управления сохраненными значениями; format или sprintf следует использовать для управления отображением значений.   -  person Carl Witthoft    schedule 21.05.2014
comment
Если у вас такие низкие p-значения, вам следует работать с их логами. Однако обычно различия между чрезвычайно низкими значениями p не являются ни значимыми, ни надежными. Таким образом, format.pval сообщает о них как "< [eps]", где eps равно .Machine$double.eps.   -  person Roland    schedule 21.05.2014
comment
Не знал о format @Andrie, спасибо за это, однако теперь кажется, что format(2.147484e-325, digits=3,scientific = T) возвращает [1] "0e+00". Правильно ли я использую format?   -  person Ashwin    schedule 21.05.2014
comment
@Roland для визуализации этих низких значений p. Я использую их журналы, но здесь я просто хотел написать их в текстовом файле в сокращенном виде.   -  person Ashwin    schedule 21.05.2014
comment
2.147484e-325 меньше, чем наименьшее положительное значение, которое может быть точно представлено в двоичном формате с плавающей запятой IEEE 754, что составляет около 4.9E-324. На результаты может повлиять небольшое количество значащих битов в представлении субнормальных чисел, чисел ниже примерно 2,22e-308.   -  person Patricia Shanahan    schedule 21.05.2014


Ответы (2)


Этот ответ основан на предположении, что числа с плавающей запятой R представлены 64-битными двоичными числами IEEE 754. Это согласуется с опубликованными результатами и по своей сути весьма вероятно.

Выполнение большого количества арифметических операций над числами с абсолютной величиной ниже примерно 2e-308 очень проблематично. Ниже этой точки точность падает с величиной. Наименьшее представимое число, около 4,9E-324, имеет в своем представлении один значащий бит. Его пара соседних чисел — 0 и примерно 1.0E-323. Любая ошибка округления уменьшит его до нуля или удвоит. Он не может столкнуться с тонкой ошибкой округления, которая влияет только на цифры с низким значением в его десятичном представлении. Точно так же round не может его немного изменить. Он может вернуть его без изменений, удвоить его, вернуть 0 или внести еще большее изменение.

См. Денормальное число для получения дополнительной информации о том, что происходит.

Решение состоит в том, чтобы по возможности избегать арифметических операций с такими крошечными числами. Один из хороших подходов, как уже отмечалось в комментариях, заключается в использовании логарифмов. Это единственный способ, если вам нужно иметь дело как с очень большими, так и с очень маленькими числами.

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

person Patricia Shanahan    schedule 21.05.2014
comment
Это правильно: см. ?.Machine и, в частности, .Machine$double.xmin информацию о R и конкретной платформе. - person Ben Bolker; 21.05.2014

Чтобы преодолеть это, вы можете использовать библиотеку Rmpfr.

Пример:

library(Rmpfr)
x <- mpfr(7.356789e-311, 80)
round(x, 313)
1 'mpfr' number of precision  1040   bits 
[1]7.36000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000008e-311
person art    schedule 01.12.2014