Logaritmo de la función Beta incompleta para grandes $\alpha,\beta$

Question

La función de R pbeta se supone que calcula la función Beta incompleta regularizada. Si la bandera log=TRUE se pasa como argumento, se devuelve el logaritmo del resultado en su lugar.

He encontrado el siguiente error numérico. Evaluando:

pbeta(0.5555555, 1925.74, 33.7179, log=TRUE)

devuelve -Inf , lo cual es obviamente erróneo. El resultado correcto es (calculado con Mathematica):

$$\log\mathrm{I}_{0.5555555}(1925.74,33.7179) = -994.767594138466967$$

Para todos los demás valores que he probado, pbeta da la respuesta correcta. He inspeccionado el código fuente de pbeta pero no puede precisar el error.

Necesito un método preciso para calcular el logaritmo de la función Beta incompleta para valores grandes de $\alpha,\beta$ (del orden de miles). ¿Puede alguien sugerir una solución para pbeta ¿o una solución? ¿O tal vez una biblioteca diferente?

Informe de errores de R: https://bugs.r-project.org/bugzilla3/show_bug.cgi?id=16332

Answer 1

La función incompleta Beta se puede escribir con la ayuda de la función hipergeométrica de Gauss: $$B_x(a,b)=\frac{1}{a}x^a{(1-x)}^b F(a+b,1,a+1,x)$$ y la FCD del $Beta(a,b)$ distribución evaluada en $x$ es $$B_x(a,b)/B(a,b).$$ Una buena implementación de la función hipergeométrica de Gauss se encuentra en el gsl paquete - una envoltura para las funciones especiales de la Biblioteca Científica Gnu. Así que usted puede escribir el log-CDF así:

library(gsl)
logpbeta <- function(x,a,b) log(hyperg_2F1(a+b,1,a+1,x)) + a*log(x)+b*log(1-x)-log(a) - lbeta(a,b)

Y da el mismo resultado que Mathematica para tu ejemplo:

> logpbeta(0.5555555, 1925.74, 33.7179)
[1] -994.7676

No sé para qué rango de parámetros el resultado es correcto. Y no he encontrado la respuesta en el manual de referencia de GNU.

Tenga en cuenta que el cero que se obtiene con pbeta parece deberse a la evaluación de $B_x(a,b)$ seguido de la transformación logarítmica:

> x <- 0.5555555 
> a <- 1925.74
> b <- 33.7179
> log(hyperg_2F1(a+b,1,a+1,x)*x^a*(1-x)^b/a)
[1] -Inf
> hyperg_2F1(a+b,1,a+1,x)
[1] 2.298761
> x^a*(1-x)^b/a
[1] 0