使用 Rmpfr 在 R 中进行高精度计算

Question

我必须计算 B(x+a, n-x+a)，其中 B( , ) 是 beta function、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 函数的关系来计算它。结果不一样。
但是先把输入转换成class"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
#------------

使用 Rmpfr 在 R 中进行高精度计算

High precision calculation in R using Rmpfr

precision

r