Работа с очень маленькими числами в R

Question

Мне нужно вычислить список очень маленьких чисел, таких как

(0,1) ^ 1000, 0,2 ^ (1200),

, а затем нормализовать их, чтобы они суммировали до одногот.е.

a1 = 0,1 ^ 1000, a2 = 0,2 ^ 1200

И я хочу вычислить a1 '= a1 / (a1 + a2), a2' = a2 (a1 + a2).

У меня проблемы с потерей памяти, так как я получаю a1 = 0.Как я могу обойти это?Теоретически я мог бы иметь дело с журналами, и тогда log (a1) = 1000 * log (0.l) был бы способом представить a1 без проблем с занижением - но для нормализации мне нужно было бы получить log (a1 + a2) -который я не могу вычислить, так как не могу представить a1 напрямую.

Я программирую на R - насколько я могу судить, нет такого типа данных, как Decimal в c #, который позволил бы вам стать лучше, чемзначение двойной точности.

Любые предложения будут оценены, спасибо

Answer 1

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

Но, чтобы сделать это в целом, вы можете использовать идею из logspace_add C-функции, которая находится под капотомR. Можно определить logxpy ( =log(x+y) ), когда lx = log(x) и ly = log(y) как:

logxpy <- function(lx,ly) max(lx,ly) + log1p(exp(-abs(lx-ly)))

Что означает, что мы можем использовать:

> la1 <- 1000*log(0.1)
> la2 <- 1200*log(0.2)

> exp(la1 - logxpy(la1,la2))
[1] 5.807714e-162

> exp(la2 - logxpy(la1,la2))
[1] 1

Эта функция может вызываться рекурсивно какхорошо, если у вас есть больше номеров.Имейте в виду, 1 все еще 1, а не 1 минус 5.807...e-162.Если вам действительно нужно больше точности, и ваша платформа поддерживает длинные двойные типы, вы можете написать все, например, на C или C ++, и вернуть результаты позже.Но если я прав, R может - на данный момент - иметь дело только с обычными удвоениями, поэтому в конечном итоге вы снова потеряете точность, когда будет показан результат.

---

РЕДАКТИРОВАТЬ:

чтобы сделать математику для вас:

log(x+y) = log(exp(lx)+exp(ly))
         = log( exp(lx) * (1 + exp(ly-lx) )
         = lx + log ( 1 + exp(ly - lx)  )

Теперь вы просто берете наибольшее значение как lx, а затем вы получаете выражение в logxpy().

РЕДАКТИРОВАТЬ 2: Зачем тогда брать максимум?Легко убедиться, что вы используете отрицательное число в exp (lx-ly).Если lx-ly становится слишком большим, тогда exp (lx-ly) возвращает Inf.Это не правильный результат.exp (ly-lx) вернет 0, что позволяет получить гораздо лучший результат:

Скажем, lx = 1 и ly = 1000, затем:

> 1+log1p(exp(1000-1))
[1] Inf
> 1000+log1p(exp(1-1000))
[1] 1000

Answer 2

Пакет Brobdingnag имеет дело с очень большими или маленькими числами, по сути, оборачивая ответ Джориса в удобную форму.

a1 <- as.brob(0.1)^1000
a2 <- as.brob(0.2)^1200
a1_dash <- a1 / (a1 + a2)
a2_dash <- a2 / (a1 + a2)
as.numeric(a1_dash)
as.numeric(a2_dash)

Answer 3

Попробуйте пакет произвольной точности mpfr.

Ryacas также может иметь произвольную точность.

Answer 4

Может быть, вы можете рассматривать a1 и a2 как дроби. В вашем примере с

a1 = (a1num/a1denom)^1000  # 1/10
a2 = (a2num/a2denom)^1200  # 1/5

вы получите

a1' = (a1num^1000 * a2denom^1200)/(a1num^1000 * a2denom^1200 + a1denom^1000 * a2num^1200)
a2' = (a1denom^1000 * a2num^1200)/(a1num^1000 * a2denom^1200 + a1denom^1000 * a2num^1200)

, который можно вычислить с помощью пакета gmp:

library(gmp)
a1 <- as.double(pow.bigz(5,1200) / (pow.bigz(5,1200)+ pow.bigz(10,1000)))