Berechnung der hypergeometrischen Funktion in R.

8

Ich habe enorme Schwierigkeiten, 2 F 1 ( a , b ; c ; z ) mit dem Paket in R zu bewerten . In meinem Fall sind die Werte von a , b , c immer positive reelle Zahlen. Trotzdem ist die hypergeometrische Funktion unglaublich empfindlich gegenüber ihren Werten. Ich suche keine extreme Präzision; Ich kann Excel verwenden, um eine grobe Schätzung der Guge-Hypergeometrie zu erhalten, die für meine Zwecke in Ordnung ist.2F1(a,b;c;z)hypergeoabc

Irgendwelche Vorschläge für eine Implementierung in R, die eine schnelle, fehlerfreie, wenn nicht sehr genaue hypergeometrische Gaußsche Berechnung positiver reeller Zahlen mit einem breiten Wertebereich ermöglicht?

Edit: Es scheint, dass es in MATLAB weit mehr Code dafür gibt als in R. Irgendwelche Gedanken darüber, warum?

benrolls
quelle
Ich hätte gedacht, ein guter Ansatz wäre es, den größten Begriff in der hypergeometrischen Summe zu finden und diesen Begriff abzuschneiden. Auf diese Weise können Sie den größten Begriff herausrechnen und mit Begriffen arbeiten, die kleiner als 1 in der Summe sind. Sie können auch die Laplace-Methode für die Integraldarstellung verwenden (geeignet parametriert, um Grenzprobleme zu vermeiden).
Wahrscheinlichkeitslogik

Antworten:

14

Wenn Sie die hypergeometrische Gauß-Funktion nicht für komplexe Werte der Parameter oder der Variablen auswerten müssen, ist es besser, das gslPaket von Robin Hankin zu verwenden .

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

Aktualisieren

Hier ist meine alternative Implementierung mit dem gmp-Paket (zumindest zum Spaß)

Stéphane Laurent
quelle
1
Vielen Dank für eine einfache, saubere Lösung für ein nicht so einfaches Problem! Es gibt nicht viel Dokumentation zu hypergeometrischen Gauß-Funktionen für R, daher ist dieser Ansatz wertvoll.
Benrolls
@benrolls Ich habe die hypergeometrische Gauß-Funktion mit dem Hypergeo-Paket, dem gsl-Paket, einer anderen von mir implementierten Methode und Mathematica erlebt. Ich garantiere, dass die oben angegebene Lösung sehr leistungsfähig ist. Ich habe bei der Verwendung nie einen Fehler gefunden.
Stéphane Laurent
ca
1

@ Stéphane Laurent's Formel oben ist großartig. Ich habe bemerkt , dass es manchmal produziert NaNs , wenn a, b, csind groß und znegativ ist - ich habe unten nicht in der Lage , die genauen Bedingungen zu fixieren. In diesen Fällen können wir eine andere hypergeometrische Transformation verwenden, die von Stéphanes alternativem Ausdruck ausgeht. Es führt zu dieser alternativen Formel:

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

Zum Beispiel:

> 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

Alle drei stimmen mit Mathematica überein Hypergeometric2F1. Diese Formel scheint für kleinere auch gut benommen a, b, c. Beachten Sie, dass es Fälle gibt, in denen diese Formel gibt NaNund Stéphane nicht . Am besten von Fall zu Fall prüfen.

pglpm
quelle