Cálculo de la función hipergeométrica en R

8

Tengo enormes dificultades para evaluar con el paquete en R. En mi caso, los valores de , , son siempre números reales positivos. Aun así, la función hipergeométrica es increíblemente sensible a sus valores. No estoy buscando precisión extrema; Puedo usar Excel para obtener una estimación aproximada de la hipergeometría de Guass que está bien para mis propósitos.2F1(a,b;c;z)hypergeoabc

¿Alguna sugerencia para una implementación en R que proporcione un cálculo hipergeométrico gaussiano rápido, libre de errores, si no súper preciso, de números reales positivos con una amplia gama de valores?

Editar: parece que hay mucho más código para esto en MATLAB que R. ¿Alguna idea de por qué?

benrolls
fuente
hubiera pensado que un buen enfoque sería encontrar el término más grande en la suma hipergeométrica y truncar alrededor de este término. De esta manera, puede factorizar el término más grande y trabajar con términos menores que 1 en la suma. También puede usar el método de Laplace en la representación integral (adecuadamente parametrizada para evitar problemas de límites).
probabilidadislogica

Respuestas:

14

A menos que necesite evaluar la función hipergeométrica de Gauss para valores complejos de los parámetros o la variable, es mejor usar el gslpaquete de Robin Hankin .

Según mi experiencia, también recomiendo evaluar solo la función hipergeométrica de Gauss para un valor de la variable que se encuentra en , y usar una fórmula de transformación para valores en .[0,1]],0]

library(gsl)
Gauss2F1 <- function(a,b,c,x){
    if(x>=0 & x<1){
        hyperg_2F1(a,b,c,x)
    }else{
            hyperg_2F1(c-a,b,c,1-1/(1-x))/(1-x)^b
        }
}

Actualizar

Aquí está mi implementación alternativa con el paquete gmp (al menos, por diversión)

Stéphane Laurent
fuente
1
¡Gracias por una solución simple y limpia a un problema no tan simple! No hay mucha documentación sobre las funciones hipergeométricas de Gauss para R, por lo que este enfoque es valioso.
benrolls
@benrolls He experimentado la función hipergeométrica de Gauss con el paquete hypergeo, el paquete gsl, otro método implementado por mí mismo y Mathematica también. Le garantizo que la solución que doy arriba es muy eficiente, nunca he encontrado un error al usarla.
Stéphane Laurent
La fórmula de Stéphane anterior es excelente, pero he notado que produce NaNs cuando es negativo con un gran valor absoluto. El uso de otra fórmula de transformación conduce a: biblioteca (gsl)C-una
pglpm
1

La fórmula anterior de @ Stéphane Laurent es excelente. Me he dado cuenta de que a veces se produce NaNs cuando a, b, cson grandes y zes negativo - No he sido capaz de precisar las condiciones precisas hacia abajo. En estos casos, podemos usar otra transformación hipergeométrica a partir de la expresión alternativa de Stéphane. Lleva a esta fórmula alternativa:

library(gsl)
Gauss2F1b <- function(a,b,c,x){
    if(x>=0 & x<1){
        hyperg_2F1(a,b,c,x)
    }else{
        hyperg_2F1(a,c-b,c,1-1/(1-x))/(1-x)^a
    }
}

Por ejemplo:

> Gauss2F1(80.2,50.1,61.3,-1)
[1] NaN
>
> Gauss2F1b(80.2,50.1,61.3,-1)
[1] 5.498597e-20
>
>
> Gauss2F1(80.2,50.1,61.3,-3)
[1] NaN
> Gauss2F1b(80.2,50.1,61.3,-3)
[1] 5.343807e-38
>
>
> Gauss2F1(80.2,50.1,61.3,-0.4)
[1] NaN
> Gauss2F1b(80.2,50.1,61.3,-0.4)
[1] 3.322785e-10

los tres están de acuerdo con Mathematica's Hypergeometric2F1. Esta fórmula parece bien comportado también para pequeños a, b, c. Tenga en cuenta que hay casos en los que esta fórmula da NaNy Stéphane no . Lo mejor es verificar caso por caso.

pglpm
fuente