R-程序设计语言中密度函数的舍入误差

oknwwptz 于 2023-03-27 发布在其他

关注(0)|答案(1)|浏览(99)

使用R，我试图计算一个值向量的似然性。其中一些值远离分布的尾部。舍入错误似乎导致结果舍入为零，使我的代码在绘制样本时抛出错误。
我尝试使用Rmpfr包来获得更高的精度，但这改变了我的变量类型。我还考虑将表达式重写为对数形式以避免乘法。然后表达式将变为：exp(log(w[X]) + dnorm(y, x[X], sigma[X], log = TRUE))
由于求幂运算，这仍然会导致函数返回零。
有没有一种方法可以用对数来公式化这个问题，允许高数值精度？我最后想要的是正确计算以下概率，总和为1。

liks <- lapply(1:3, function(X) {
    w[X] * dnorm(y, x[X], sigma[X])
  })

  probs1 <- liks[[1]]/(liks[[1]] + liks[[2]] + liks[[3]])
  
  probs2 <- liks[[2]]/(liks[[1]] + liks[[2]] + liks[[3]])
  
  probs3 <- liks[[3]]/(liks[[1]] + liks[[2]] + liks[[3]])

编辑：增加日志代码数值示例

下面是一些参数和示例值，它们开始抛出日志代码错误。它们不会抛出NAN。相反，三个条目中的每一个都是零：

$mu
[1] 0.7323412910 0.7742235621 0.4863889347

$w
[1] 0.008464 0.083536 0.908000

$sigma
[1] 0.08209500030 0.08166088502 0.09168991045

Observation values:
c(4.667935371,  5.654500961,  4.383309364,  4.396201611,  4.452524185,  4.441100597,  4.890487194,  4.416962624,  5.241273880,  4.347382069,  4.867616177,  4.895996094,  4.592288494, -3.612523079,  4.817468166,  4.783963203,  4.541391850,  4.709537983,  5.227987289,  5.585811138,  4.497674942,  4.989979267,  4.489729881)

在我的数据集中，所有的观察值都是极端值/异常值。这可以解释为什么它们被分配的概率如此之小，以至于它们被四舍五入为零。

r

来源：https://stackoverflow.com/questions/75842424/rounding-error-in-density-functions-of-r-programming-language

1条答案

按热度按时间

kt06eoxx1#

你应该在整个过程中使用日志，不要这么快就使用指数。例如，这里有一个变化：

logliks <- lapply(1:3, function(X) {
  log(w[X]) +  dnorm(y, x[X], sigma[X], log = TRUE)
 })

现在，要计算像probs1这样的表达式，您需要将分子和分母除以liks值中的最大值，即计算

(liks[[1]]/biggest)/(liks[[1]]/biggest + liks[[2]]/biggest + liks[[3]]/biggest)

但都是对数尺度的

logbiggest <- max(as.numeric(logliks))
logprobs1 <- (logliks[[1]] - logbiggest) - 
  log( exp( logliks[[1]] - logbiggest ) 
    + exp(logliks[[2]] - logbiggest) 
    + exp(logliks[[3]] - logbiggest) )

对于logprobs2和logprobs3也是如此。由于logbiggest等于logliks中的一个，因此这些指数中的一个将等于1.0，然后其他指数是否下溢无关紧要。
编辑以添加：数值例子
您在问题中添加了数据。下面是使用您的数据进行的完整计算。我没有得到任何零概率，但大多数概率都非常小：

x <- c(0.7323412910, 0.7742235621, 0.4863889347)
w <- c(0.008464, 0.083536, 0.908000)
sigma <- c(0.08209500030, 0.08166088502, 0.09168991045)

y <- c(4.667935371,  5.654500961,  4.383309364,  4.396201611,  4.452524185,  4.441100597,  4.890487194,  4.416962624,  5.241273880,  4.347382069,  4.867616177,  4.895996094,  4.592288494, -3.612523079,  4.817468166,  4.783963203,  4.541391850,  4.709537983,  5.227987289,  5.585811138,  4.497674942,  4.989979267,  4.489729881)

probs <- matrix(NA, length(y), 3)

for (i in seq_along(y)) {
  logliks <- lapply(1:3, function(X) {
    log(w[X]) +  dnorm(y[i], x[X], sigma[X], log = TRUE)
  })
  
  logbiggest <- max(as.numeric(logliks))
  
  logprobs1 <- (logliks[[1]] - logbiggest) - 
    log( exp( logliks[[1]] - logbiggest ) 
         + exp(logliks[[2]] - logbiggest) 
         + exp(logliks[[3]] - logbiggest) )
  
  logprobs2 <- (logliks[[2]] - logbiggest) - 
    log( exp( logliks[[1]] - logbiggest ) 
         + exp(logliks[[2]] - logbiggest) 
         + exp(logliks[[3]] - logbiggest) )
  
  logprobs3 <- (logliks[[3]] - logbiggest) - 
    log( exp( logliks[[1]] - logbiggest ) 
         + exp(logliks[[2]] - logbiggest) 
         + exp(logliks[[3]] - logbiggest) )
  
  probs[i, ] <- exp(c(logprobs1, logprobs2, logprobs3))
}

head(probs)
#>              [,1]         [,2] [,3]
#> [1,] 4.006975e-50 9.059815e-44    1
#> [2,] 1.964023e-93 2.172199e-87    1
#> [3,] 6.103586e-40 1.274113e-33    1
#> [4,] 2.221979e-40 4.668064e-34    1
#> [5,] 2.538854e-42 5.467829e-36    1
#> [6,] 6.336084e-42 1.358276e-35    1

创建于2023-03-26带有reprex v2.0.2
由于四舍五入，最后一个等于1;这是有道理的，因为y的值比均值大得多，第三个分布的方差最大：所以模型预测所有这些异常值都可能来自该分布。将sigma乘以10，你会得到更少的极端概率。

赞(0）回复(0）举报 2023-03-27

我来回答

R-程序设计语言中密度函数的舍入误差

1条答案

相关问题

热门标签

最新问答