The beta function (not the beta distribution) beta(a,b) can be computed when a or b or both are negative. However, in the implementation of R, only non-negative arguments are allowed.
Is there a way around this?
A common definition of beta(a,b) is gamma(a)*gamma(b)/gamma(a+b); hence, we could use lgamma() to avoid overflow:
exp(lgamma(a)+lgamma(b)-lgamma(a+b))
but for some a or some b, say, when a <- -2.5; b <- 3, the solution to lgamma(a) should be a complex number yet R only returns the real part...
pracma::gammaz or gsl::lngamma_complex seem to work (and agree with Wolfram Alpha).
## convert values to complex (not sure if this is necessary)
a <- -2.5+0i; b <- 3+0i
library(pracma)
Re(gammaz(a)*gammaz(b)/(gammaz(a+b))) ## -1.0667
Better (since as you point out it avoids overflow problems when a and b are large):
library(gsl)
exp((lngamma_complex(a)+lngamma_complex(b))-lngamma_complex(a+b))
## [1] -1.066667-1.30629e-16i
It's not obvious to me when the Beta function is real-valued; if you were confident that the result will be real, you can use Re() to drop the imaginary part of the result. Otherwise this answer might come in handy for dropping floating-point fuzz in the imaginary part ...