Обработка очень малых чисел в соотношении и как сохранить экспоненциальную стоимость

Question

В настоящее время я использую R версии 3.4.4 (2018-03-15) с R Studio.

Мне нужно рассчитать соотношение двух значений.И у меня есть проблемы с некоторыми случаями:

• Числитель может иметь очень малое значение: exp (-2408.9), что R приближается к 0.
• Знаменатель также: exp (-2405) рассчитывается как 0 - это R.

Когда вычисляется отношение, я получаю NaN (из-за 0/0).

Первое решение :

Я использую библиотеку Brobdingnag , которая позволяет сохранить число в качестве экспоненты, и, наконец, получим, что соотношение на самом деле: exp (-3.8987) = 0.02026725

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

Вопросы для другого решения:

Вы слышали о другой библиотеке, которая работает с очень маленькими (или большими) значениями?

Я бы хотел сохранить мой числитель и знаменатель в экспоненциальном выражении, пока не будет выполнено деление, но я понятия не имею, как это сделать.Конечно, мой числитель и знаменатель являются векторами, которые я делю, как только они оба вычисляются.(Я не могу получить знаменатель без вектора числителя) Есть ли способ «заставить» R сохранить значение как exp вместо целого числа (и 0 ...)?

Заранее спасибо за любую помощь.

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

Вот соотношение, которое я должен рассчитать:

https://ibb.co/dFHx4z

Я не уверен, что могу использовать трюк: exp(x) / exp (y) = exp (xy), потому что у меня есть сумма в деноме.Вот почему мне нужна формула exp, пока я не сделаю соотношение ... Значения внутри exp - очень большое отрицательное число, и exp из этих чисел равняется 0. Кроме того, я попытался преобразовать числитель как log, чтобы я мог получить журнал первыхПрошлая + вторая часть (без опыта), но иногда первая часть числителя (1 / sqrt ...) имеет малый размер, и журнал ее возвращает Inf ..

Я думаю, что есть способ, но я могу' Найди это.

Спасибо за все ответы, кстати!

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

####### Fonction that calculate the density (with brobdingnag package) :

density <- function(nc,yc,X,beta,sig,k){

    # n_c is a vector of integer 
    # y_c is a vector of numeric 
    # X is a matrix 
    # beta is a vector of numeric 
    # sigma is a value

res<-as.brob((1/(2*pi*sig[k])))^(nc/2)*exp(as.brob(-(1/(2*sig[k]))*t(yc-(X %*% beta[,k])) %*% (yc-(X %*% beta[,k]))))
return(res)
}

####### Code for calculation of the ratio :

# n_c[c] : num [1] 340
# y_c[c] : num [1:340] 1.279 0.777 1.069 0.864 1.56 ...
# X[c] : num [1:340, 1:11] 1 1 1 1 1 1 1 1 1 1 ... (matrix of 0 and 1)
# beta : num [1:11, 1:2] 1.542 -0.226 -0.145 -0.438 -0.201 ...
# sigma : num [1:2] 21.694381  4.267277
# lambda : num [1] 0.5

# Numerator :

num_tau<-sapply(1:100,function(c){
        sapply(1:4,function(k){
            lambda[k]*density(n_c[c], y_c[c],X[c],beta,sigma,k)
        })
    })

# Denominator :

denom_tau<-list()
for (c in 1:100){
    val<-0
    for (k in 1:4){
        val<-val+num_tau[k,c][[1]]
    }
denom_tau[[c]]<-val
}

# Ratio :
for (l in 1:4){
    for (c in 1:100){
        tau[l,c]<-as.numeric(num_tau[l,c][[1]]/denom_tau[[c]])
    }
}

Answer 1

Как рекомендует @minem, вы можете использовать пакет Rmpfr .Вот один из способов применить его к вашему случаю.

Сначала переместите множители внутри экспоненты числителя, используя тот факт, что a * exp (b) = exp (b + log (a)).Затем переписайте функцию density, чтобы вычислить числитель журнала:

log_numerator <- function(nc, yc, X, beta, sig, k, lambda){
  v <- yc - X %*% beta[,k]
  res <- -sum(v*v)/(2*sig[k]) - (nc/2)*log(2*pi*sig[k]) + log(lambda[k])
  drop(res)
}

Обратите внимание, что lambda теперь передается этой функции.Также обратите внимание, что мы можем вычислить точечное произведение вектора Y - X * бета, как показано.

Теперь мы можем сгенерировать некоторые данные.Здесь я исправляю c и просто k = 1: 2.

set.seed(1)
n_c <- 340
y_c <- rnorm(340)
dat <- data.frame(fac = sample(letters[1:11], 340, replace = TRUE)
X_c <- model.matrix(~ fac, data = dat)
beta <- matrix(runif(22, -10, 10), 11, 2)
sigma <- c(21.694381,  4.267277)
lambda <- c(0.5, 0.5)

Используя вашу функцию плотности, мы имеем

x1 <- lambda[1] *density(n_c, y_c,X_c,beta,sigma,1)
y1 <- lambda[2] *density(n_c, y_c,X_c,beta,sigma,2)
x1
# [1] +exp(-1738.4)
y1
# [1] +exp(-1838.7)
as.numeric(y1/sum(x1, y1))
# [1] 2.780805e-44

Используя функцию числового логарифма, мы имеем

p <- 40
x <- mpfr(log_numerator(n_c, y_c,X_c,beta,sigma,1, lambda), p)
y <- mpfr(log_numerator(n_c, y_c,X_c,beta,sigma,2, lambda), p)
x
# 1 'mpfr' number of precision  40   bits 
# [1] -1738.379327798
y
# 1 'mpfr' number of precision  40   bits 
# [1] -1838.67033143
exp(y)/sum(exp(x), exp(y))
# 1 'mpfr' number of precision  53   bits 
# [1] 2.780805017186589e-44

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

Вы также можете повысить эффективность, используя больше векторизации.Например, мы можем векторизовать log_numerator по k:

log_numerator2 <- function(nc, yc, X, beta, sig, lambda){
  M <- yc - X %*% beta
  res <- -colSums(M*M)/(2*sig) - (nc/2)*log(2*pi*sig) + log(lambda)
  drop(res)
}
z <- log_numerator2(n_c, y_c, X_c, beta, sigma, lambda)
z
# [1] -1738.379 -1838.670

Теперь предположим, что у нас есть числовые логарифмы в матрице ac на k, для иллюстрации предположим, что все c имеют те же значения, что и z,

log_num <- mpfr(matrix(z, byrow = TRUE, 3, 2), p)

Вы можете вычислить отношения следующим образом

num <- exp(log_num)
denom <- apply(num, 1, sum) # rowSums not implemented for mpfr
num/denom
# 'mpfrMatrix' of dim(.) =  (3, 2) of precision  53   bits 
#     [,1]              [,2]                 
# [1,] 1.000000000000000 2.780805017186589e-44
# [2,] 1.000000000000000 2.780805017186589e-44
# [3,] 1.000000000000000 2.780805017186589e-44

Answer 2

Если оба значения нужно предварительно экспонировать, вы можете использовать формулу:

e ^ x / e ^ y = e ^ (xy)

В противном случае вы можете попробовать Rmpfr пакет.

Пример:

require(Rmpfr)
p = 40
x <- mpfr(-2408.9, p)
y <- mpfr(-2405, p)
exp(x)/exp(y)
# 1 'mpfr' number of precision  40   bits 
# [1] 0.02024191147598