处理非常小的比率数字以及如何保持指数值

Question

我目前在 R Studio 中使用 R 版本 3.4.4 (2018-03-15)。

我需要计算两个值的比率。我在某些情况下遇到问题：

分子可以是非常小的值：exp(-2408.9)，R近似于0。
分母也 : exp(-2405) 计算为 0 是 R.

计算比率时，我得到一个 NaN（因为 0/0）。

第一个解 :

我使用 Brobdingnag 库，它允许将数字保持为指数，最后得到的比率实际上是：exp(-3.8987) = 0.02026725

但是，使用库 profvis 检查我的代码的性能时，我发现尽管 Brobdingnag 库对我来说非常有用，但它让我付出了代价很多在性能方面。而且我不能保留这个解决方案，因为我必须对我的算法进行大量模拟。

其他解决方案的问题：

您是否听说过其他库可以处理非常小（或大）的值？

我想在进行除法之前将我的分子和分母保留在指数表达式中，但我不知道该怎么做。因为当然，我的分子和分母是向量，一旦它们都被计算出来，我就将它们相除。（没有分子向量我无法获得分母）有没有办法 "force" R 将值保留为 exp 而不是整数（和 0...）？

提前感谢您的帮助。

编辑：

这是我必须计算的比率：

https://ibb.co/dFHx4z

我不确定我是否可以使用这个技巧：exp(x)/exp(y) = exp(x-y) 因为我在 denom 中有一个总和。这就是为什么在计算比率之前我需要 exp 公式的原因...... exp 中的值是非常大的负数，这些数字的 exp 为 0。另外，我尝试将分子转换为对数，因此我可以得到第一部分 + 第二部分（没有 exp）的对数，但有时，分子的第一部分(1/sqrt...) 太小了，记录下来 returns 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

如果两个值之前都需要取幂，那么可以使用公式：

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

否则你可以试试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

Answer 2

正如@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*beta 的点积，如图所示。

现在我们可以生成一些数据了。在这里我修复了 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 可以用来产生相同的结果，但是如果没有更好的测试代码就很难检查时间。

您还可以通过使用更多矢量化来提高效率。例如。我们可以在 k 上对 log_numerator 进行向量化：

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

现在假设我们在 c x 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

处理非常小的比率数字以及如何保持指数值

Handling very small numbers in ratio and how to keep exponential value

performance

r

exponential