r-Rcpp中对数CDF值的欠流处理

r - Treatment of Underflow from Log CDF Values in Rcpp

本文关键字：处理 CDF r-Rcpp      更新时间：2023-10-16

我目前正在使用Rcpp编写一个R包，需要评估使用beta分布的日志可能性。根据数据的大小（影响贝塔中的a和b参数），一些左尾贝塔分布概率最终非常小。由于我正在获取日志，我最终会遇到下溢问题。例如：

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")

如果你运行这个代码，你会得到警告：

警告消息：在函数（）中：pbeta（*，log.p=TRUE）->bpser（a=6403，b=21，x=0.8，…）下溢到-Inf

我知道我可能可以在函数本身中使用if语句将-Inf替换为-DBL_MAX来解决后续的计算问题，但这可能并不能消除警告消息。你知道如何抑制警告消息（或者更优雅地处理这个问题）吗？提前谢谢。

我们应该提供与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（当然）应该尽可能精确——这就是我们从像马丁·梅克勒这样痴迷细节的严肃统计学家（TM）的辛勤工作中得到的。我认为您应该在r-devel上尝试这个问题，因为我没有看到Rcpp在这里做错了什么——它只是将值封送至r本身使用的函数。

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

相关文章：