Rmpfr在R中的高精度计算

Question

我必须计算B(x+a, n-x+a)，其中B( , )是贝塔函数，0 < a < .5，n = 10^8和x = 10^7。

R只是喷出0，这破坏了我的计算。我尝试使用Rmpfr包，但是在beta函数中使用mpfr对象时出现了一个错误，如下所示：

beta(Rmpfr::mpfr(.3, 32), Rmpfr::mpfr(.4, 32))

Error in beta(Rmpfr::mpfr(0.3, 32), Rmpfr::mpfr(0.4, 32)) : non-numeric argument to mathematical function

在Rmpfr包中是否有使用此函数或扩展的方法？

Answer 1

这里有个办法。

定义一个计算B(x + a, n - x + a)的函数，或者使用lbeta或与Gamma函数的关系来计算它。结果是不同的。

但是首先要转换"mpfr"类中的输入。

library(Rmpfr)

fun1 <- function(x, a, n, precBits = 128){
  x <- mpfr(x, precBits = precBits)
  y <- lbeta(x + a, n - x + a)
  exp(y)
}
fun2 <- function(x, a, n, precBits = 128){
  x <- mpfr(x, precBits = precBits)
  y <- lgamma(x + a) + lgamma(n - x + a) - lgamma(n + 2*a)
  exp(y)
}

x <- 10^7
n <- 10^8

a <- seq(0, 0.5, by = 0.1)
a[1] <- a[1] + .Machine$double.eps

y1 <- sapply(a, function(.a) fun1(x, .a, n))
y2 <- sapply(a, function(.a) fun2(x, .a, n))

identical(y1, y2)    # FALSE

for(i in seq_along(y1)){
  print(y1[[i]])
  print(y2[[i]])
  cat("------------\n")
}
#'mpfr1' 5.908917437507173802740403605476917652884e-14118178
#'mpfr1' 5.908916502709463876883575806175392305756e-14118178
#------------
#'mpfr1' 4.644427318909735435193584822408613288295e-14118178
#'mpfr1' 4.644427671609397947912541833050543093524e-14118178
#------------
#'mpfr1' 3.650534190755993708699375955051812478253e-14118178
#'mpfr1' 3.650534453825956612660829361117675142023e-14118178
#------------
#'mpfr1' 2.869331130039816933725658480797135290929e-14118178
#'mpfr1' 2.869331326837006304711880791004313346298e-14118178
#------------
#'mpfr1' 2.255303117148781009056940914963421818211e-14118178
#'mpfr1' 2.255303264892447200718169160444143102796e-14118178
#------------
#'mpfr1' 1.772675206631663936649823041615115643639e-14118178
#'mpfr1' 1.772675318013218204950148512940582629174e-14118178
#------------