Computation of hypergeometric function in R

Question

I'm having tremendous difficulty evaluating $_2F_1(a,b;c;z)$ with the hypergeo package in R. In my case, values of $a$, $b$, $c$ are always positive real numbers. Even so, the hypergeometric function is incredibly sensitive to their values. I am not looking for extreme precision; I can use Excel to get a rough estimation of the Guass hypergeometric that is fine for my purposes.

Any suggestions for an implementation in R that will give a fast, error-free, if not super accurate Gaussian hypergeometric computation of positive real numbers with a wide range of values?

Edit: it seems there is far more code for this in MATLAB than R. Any thoughts as to why?

i would have thought a good approach would be to find the largest term in the hypergeometric sum and truncate around this term. this way you can factor out the largest term and work with terms less than 1 in the sum. you could also use laplace method on the integral representation (suitably parameterised to avoid boundary issues). — probabilityislogic, Aug 01 '12 at 02:16

Stéphane Laurent · Accepted Answer · 2013-03-08T18:58:33.420

14

Unless you need to evaluate the Gauss hypergeometric function for complex values of the parameters or the variable, it is better to use Robin Hankin's gsl package.

Based on my experience I also recommend to only evaluate the Gauss hypergeometric function for a value of the variable lying in $[0,1]$, and use a transformation formula for values in $]-\infty, 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
        }
}

Update

Here is my alternative implementation with the gmp package (at least, for fun)

edited Mar 08 '13 at 18:58

answered Aug 01 '12 at 05:03

Stéphane Laurent

17,425
5
59
101

1

Thanks for a simple, clean solution to a not so simple problem! There isn't a lot of documentation on gauss hypergeometric functions for R so this approach is valuable. – benrolls Aug 01 '12 at 13:39
@benrolls I have experienced the Gauss hypergeometric function with the hypergeo package, the gsl package, another method implemented by myself, and Mathematica too. I guarantee that the solution I give above is very performant, I have never find a bug when using it. – Stéphane Laurent Aug 01 '12 at 14:31
Stéphane's formula above is great, but I've noticed that it produces NaNs when $c-a$ is negative with large absolute value. Use of another transformation formula leads to: library(gsl) – pglpm Aug 06 '18 at 17:08

score 1 · Answer 2 · answered Aug 06 '18 at 19:04

@Stéphane Laurent's's formula above is great. I've noticed that it sometimes produces NaNs when a, b, c are large and z is negative – I haven't been able to pin the precise conditions down. In these cases we can use another hypergeometric transformation starting from Stéphane's alternative expression. It leads to this alternative formula:

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

For example:

> 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

all three agreeing with Mathematica's Hypergeometric2F1. This formula seems well behaved also for smaller a, b, c. Note that there are cases in which this formula gives NaN and Stéphane's doesn't, though. Best to check case by case.

Computation of hypergeometric function in R

2 Answers2

Update