Calcul de la fonction hypergéométrique dans R

8

J'ai énormément de difficulté à évaluer avec le package en R. Dans mon cas, les valeurs de , , sont toujours des nombres réels positifs. Même ainsi, la fonction hypergéométrique est incroyablement sensible à leurs valeurs. Je ne recherche pas une précision extrême; Je peux utiliser Excel pour obtenir une estimation approximative de l'hypergéométrique de Guass qui convient à mes besoins.2F1(a,b;c;z)hypergeoabc

Avez-vous des suggestions pour une implémentation en R qui donnerait un calcul hypergéométrique gaussien rapide, sans erreur, sinon super précis, de nombres réels positifs avec une large gamme de valeurs?

Edit: il semble qu'il y ait beaucoup plus de code pour cela dans MATLAB que R. Des réflexions sur pourquoi?

r hypergeometric enroule
la source
j'aurais pensé qu'une bonne approche serait de trouver le plus grand terme dans la somme hypergéométrique et de le tronquer autour de ce terme. de cette façon, vous pouvez factoriser le terme le plus grand et travailler avec des termes inférieurs à 1 dans la somme. vous pouvez également utiliser la méthode laplace sur la représentation intégrale (convenablement paramétrée pour éviter les problèmes de limites).
Probabilislogic

Réponses:

14

À moins que vous n'ayez besoin d'évaluer la fonction hypergéométrique de Gauss pour les valeurs complexes des paramètres ou de la variable, il est préférable d'utiliser le gslpackage de Robin Hankin .

Sur la base de mon expérience, je recommande également d'évaluer uniquement la fonction hypergéométrique de Gauss pour une valeur de la variable située dans et d'utiliser une formule de transformation pour les valeurs dans .[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
        }
}

Mise à jour

Voici mon implémentation alternative avec le paquet gmp (au moins, pour le plaisir)

Stéphane Laurent
la source
1
Merci pour une solution simple et propre à un problème pas si simple! Il n'y a pas beaucoup de documentation sur les fonctions hypergéométriques de Gauss pour R, cette approche est donc précieuse.
s'inscrit
@benrolls J'ai expérimenté la fonction hypergéométrique de Gauss avec le paquet hypergeo, le paquet gsl, une autre méthode implémentée par moi-même, et Mathematica aussi. Je vous garantis que la solution que je donne ci-dessus est très performante, je n'ai jamais trouvé de bug lors de son utilisation.
Stéphane Laurent
La formule de Stéphane ci-dessus est excellente, mais j'ai remarqué qu'elle produit des NaN lorsque est négatif avec une grande valeur absolue. L'utilisation d'une autre formule de transformation conduit à: bibliothèque (gsl)c−a
pglpm
1

@ La formule de Stéphane Laurent ci-dessus est excellente. Je l' ai remarqué qu'il produit parfois NaNs quand a, b, csont grandes et zest négative - je ne l' ai pas été en mesure de cerner les conditions précises vers le bas. Dans ces cas, nous pouvons utiliser une autre transformation hypergéométrique à partir de l'expression alternative de Stéphane. Cela conduit à cette formule alternative:

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
    }
}

Par exemple:

> 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

tous trois d'accord avec Mathematica Hypergeometric2F1. Cette formule semble bien se comportait aussi pour les plus petits a, b, c. Notez qu'il existe des cas où cette formule donne NaNet Stéphane de ne pas , cependant. Mieux vaut vérifier au cas par cas.

pglpm
la source