在 R 中处理非常精确的数字-解网

问：

我正在使用 R 尝试求解公式，其中和是标准正态分布函数。我遇到的问题是 g 的某些值非常小或在 R 中编码为 0，因此给定的解决方案是 Inf。下面的代码行给出了 g 在应该为正数时被编码为 0 的示例。

pnorm(29.1,0,1) - pnorm(29,0,1) #this returns a value of 0

鉴于这个问题，我有两个问题。首先，有没有一个数学技巧可以用来从上面的代码中得到一个正数？我知道有一个参数返回概率的日志。也许这会派上用场？pnormlog.p=TRUE

我的第二个问题是，如果我有一个正但非常小的 g，是否有数学技巧可以确保 f 和 g 之比的对数是有限的？

R 数学精度

密度函数（在 R 中）是 CDF （）的导数。这意味着当接近 0.因此，对于小，我们可以近似.因此（大约）和 .dnorm()pnorm()dnorm(x) = (pnorm(x + delta) - pnorm(x)) / deltadeltadeltadnorm(x) * delta = pnorm(x + delta) - pnorm(x)g(x) = dnorm(x) * dlog(g(x)) = log(dnorm(x)) + log(d)

在计算中使用为 .dnorm(29, log = TRUE) + log(0.1)log(g)

2赞 Ben Bolker 10/6/2022 #2

正如在回答和评论中指出的那样，您实际上只需要计算 .log(g)

一种解决方案是使用对数空间加法/减法，它使用一个很好的技巧来准确计算。log(exp(x) ± exp(y))

DPQ::logspace.sub(pnorm(29.1, log.p = TRUE), pnorm(29, log.p = TRUE))
## [1] -424.8435

@Alex指出，如果我们通过反转符号来观察（对称的）相反的尾巴，以便我们从另一个非常小的数字中减去一个非常小的数字，那么四舍五入误差就不是问题：

log(pnorm(-29,0,1) - pnorm(-29.1,0,1))
##  -424.8435

比较@GregorThomas的基于导数的估计值：

dnorm(29, log = TRUE) + log(0.1)
## [1] -423.7215

使用软件包检查这些解决方案：Rmpfr

library(Rmpfr)
x1 <- mpfr(29.1, 1000)  ### 1000-bit precision
x2 <- mpfr(29.0, 1000)
log(pnorm(x1) - pnorm(x2))
## 1 'mpfr' number of precision  1000   bits 
## [1] -424.843525917118620331088762368642820931865149932737039093516330881321230653835360169156688142607902967369067320263733953321675429119215305576014352760866238194504221709021981332048618426550583180796866417925580695487036011142823944123174363138216647098415472494704193913895461975221571930136768807682871

所以看起来@GregorThomas的解决方案非常接近，但其他两个解决方案更好。

在 R 中处理非常精确的数字

Handling very precise numbers in R

评论

评论