r - Rcpp 中对数 CDF 值下溢的处理

Question

我目前正在使用 Rcpp 编写 R 程序包，并且需要评估使用 beta 分布的对数似然。根据数据的大小(影响 beta 中的 a 和 b 参数)，一些左尾 beta 分布概率最终非常小。由于我正在记录日志，因此我遇到了下溢问题。例如:

library(Rcpp)
library(inline)
src <- "double x = Rf_pbeta(0.8, (double) 6403, (double) 21, 1, 1);
        return(Rcpp::wrap(x));"
fun <- cxxfunction(signature(), src, plugin = "Rcpp")

如果您运行这段代码，您会收到警告:

Warning message: In fun() : pbeta(*, log.p=TRUE) -> bpser(a=6403, b=21, x=0.8,...) underflow to -Inf

我知道我可能可以使用 if 语句将函数本身中的 -Inf 替换为 -DBL_MAX 以修复后续计算问题，但大概这不会消除警告消息。关于如何抑制警告消息(或更优雅地处理此问题)的任何想法？提前致谢。

Answer 1

我们应该提供与 R 本身完全相同的界面。由于您已经在使用它，我们可以在 R 中完成它:

R> pbeta(seq(0.2, 0.8, by=0.1), 6403, 21, log=TRUE)
[1] -10176.71  -7583.18  -5744.24  -4319.09  -3156.15  -2174.90      -Inf
Warning message:
In pbeta(seq(0.2, 0.8, by = 0.1), 6403, 21, log = TRUE) :
  pbeta(*, log.p=TRUE) -> bpser(a=6403, b=21, x=0.8,...) underflow to -Inf
R> 

所以 R 应该(当然)在这里尽可能精确——这就是我们从像 Martin Maechler 这样痴迷于细节的严肃统计学家 (TM) 的辛勤工作中得到的结果。我认为你应该在 r-devel 上尝试这个问题，因为我没有看到 Rcpp 在这里做错了什么——它只是将非常值编码到 R 本身使用的非常函数。

最后，您使用了 Rf_pbeta，这是我最不喜欢的习惯用法。考虑标量 R::pbeta 或(基于 Rcpp 糖的)矢量化 Rcpp::pbeta。