High precision calculation in R using Rmpfr

Question

I have to calculate B(x+a, n-x+a) where B( , ) is the beta function, 0 < a < .5, n = 10^8 and x = 10^7.

R just spits out 0 and this ruins my calculations further down the line. I tried using the Rmpfr package, but I get an error when using an mpfr object in the beta function like:

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

Is there a way to use this function or an extension in the Rmpfr package?

Answer 1

Here is a way.
Define a function to compute B(x + a, n - x + a) and either use lbeta or the relation with the Gamma function to compute it. The results are different.
But first transform the input in 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
#------------