This paper introduces the hypergeo package of R routines for numerical calculation of hypergeometric functions. The package is focussed on efficient and accurate evaluation of the Gauss hypergeometric function over the whole of the complex plane within the constraints of fixed-precision arithmetic. The hypergeometric series is convergent only within the unit circle, so analytic continuation must be used to define the function outside the unit circle. This short document outlines the numerical and conceptual methods used in the package; and justifies the package philosophy, which is to maintain transparent and verifiable links between the software and Abramowitz and Stegun (1965). Most of the package functionality is accessed via the single function hypergeo()
, which dispatches to one of several methods depending on the value of its arguments. The package is demonstrated in the context of game theory.
The geometric series
where
where
where it is understood that
is implemented in the package as genhypergeo_series()
and operates by
repeatedly incrementing the upper and lower index
vectors
In most cases of practical interest one finds that
Michel and Stoitsov (2008) state that physical applications are “plethora”; examples
would include atomic collisions (Alder et al. 1956), cosmology (Cruz-Dombriz and Dobado 2006),
and analysis of Feynman diagrams (Davydychev and Kalmykov 2004). In addition,
naturally-occuring combinatorial series frequently have a sum
expressible in terms of hypergeometric functions (Petkovšek et al. 1997). One
meets higher-order hypergeometric functions occasionally; the
hypergeometric distribution, for example, has a cumulative distribution
function involving the
There are two other numerical implementations for the hypergeometric
function for R: the gsl
package (Hankin 2006b), a wrapper for the Gnu Scientific Library,
although this does not cover complex values (Galassi et al. 2013); and the
appell
package (Bove et al. 2013) which implements the Gauss hypergeometric function
as hyp2f1()
.
Outside the R world, there are several proprietary implementations but the evaluation methodology is not available for inspection. Open-source implementations include that of Sage (Stein et al. 2015) and (Maxima 2014). The hypergeo package is offered as an R-centric suite of functionality with an emphasis on multiple evaluation methodologies, and transparent coding with nomenclature and structure following that of Abramowitz and Stegun (1965). An example is given below in which the positions of the cut lines may be modified.
The hypergeometric function’s series representation, namely
has unit radius of convergence by the ratio test [NB: equations with three-part numbers, as (2) above, are named for their reference in Abramowitz and Stegun (1965)]. However, the integral form
due to Gauss, furnishes analytic continuation; it is usual to follow
Riemann and define a cut along the positive real axis from f15.3.1()
in the package and exhibits surprisingly
accurate evaluation.
Gauss also provided a continued fraction form for the hypergeometric
function (implemented as hypergeo_contfrac()
in the package) which has
superior convergence rates for parts of the complex plane at the expense
of more complicated convergence properties (Cuyt et al. 2008).
The hypergeo package provides some functionality for the
hypergeometric function. the emphasis is on fast vectorized R-centric
code, complex
The majority of the package functionality is accessed via the
hypergeo()
function whose behaviour is discussed below.
Observing the slow convergence of the series representation (2), the complex behaviour of the continued fraction representation, and the heavy computational expense of the integral representation (3), it is clear that non-trivial numerical techniques are required for a production package.
The package implements a generalization of the method of Forrey (1997) to
the complex case. It utilizes the observation that the ratio of
successive terms approaches
The primary argument in equations (4)–(9) is a
member of the set
genhypergeo_series()
; see Figure 1 for a diagram
showing which parts of the complex plane use which transformation.
Given the appropriate transformation, the right hand side is evaluated
using direct summation. If genhypergeo_series()
which evaluates the generalized hypergeometric
function, Equation (1); the R implementation
uses multiplication by repeatedly incremented upper and lower
indices
Thus for example if f15.3.8()
:
> require("hypergeo")
> f15.3.8
function(A, B, C, z, tol = 0, maxiter = 2000) {
jj <- i15.3.8(A, B, C)
jj[1] * (1-z)^(-A) * genhypergeo(U = c(A, C-B), L = A-B+1, z = 1/(1-z), tol = tol,
maxiter = maxiter) + jj[2] * (1-z)^(-B) * genhypergeo(U = c(B, C-A), L = B-A+1,
z = 1/(1-z), tol = tol, maxiter = maxiter)
}
(slightly edited in the interests of visual clarity). This is a typical
internal function of the package and like all similar functions is named
for its equation number in (Abramowitz and Stegun 1965). Note the helper function
i15.3.9()
, which calculates the Gamma coefficients of the two
hypergeometric terms in the identity. This structure allows transparent
checking of the code.
The hypergeometric differential equation
together with a known value of
The package includes functionality for solving
equation (10) using ode()
from the deSolve
package (Soetaert et al. 2010):
> f15.5.1(
+ A = 1.1, B = 2.2, C = 3.5, z = 3+1i, startz = 0.5i,
+ u = function(u) straight(u, 0.5i, 3+1i),
+ udash = function(u) straightdash(u, 0.5i, 3+1i))
[1] -0.5354302+0.7081344i
> f15.5.1(
+ A = 1.1, B = 2.2, C = 3.5, z = 3+1i, startz = 0.5i,
+ u = function(u) semicircle(u, 0.5i, 3+1i, FALSE),
+ udash = function(u) semidash(u, 0.5i, 3+1i, FALSE))
[1] -1.395698-0.043599i
> hypergeo(1.1, 2.2, 3.5, 3+1i)
[1] -0.5354302+0.7081338i
See how the different integration paths give different results; the
straight path value matches that of hypergeo()
. The package also
provides hypergeo_press()
, which is somewhat more user-friendly but
less flexible, and uses the method recommended by Press et al. (1992).
The series methods detailed above are not applicable for all values of
the parameters hypergeo_cover1()
which uses (4) through (9) but
with (6) replaced with suitable limiting forms such as
This equation is comparable to (6) in terms of computational
complexity but requires evaluation of the digamma function genhypergeo_series()
but includes a runtime option
which specifies whether to evaluate
A similar methodology is used for the case hypergeo_cover2()
.
However, the case hypergeo_cover3()
which uses formulae
taken from the Wolfram functions site (Wolfram 2014). For example
w07.23.06.0026.01()
gives a straightforwardly implementable numerical
expression for
In all these cases, the limiting behaviour is problematic. For example,
consider a case where f15.3.1()
or the continued fraction
form genhypergeo_contfrac()
.
All the above methods fail when
These points were considered by (Buhring 1987) who presented a
computational method for these values; however, his method is not
suitable for finite-precision arithmetic (a brief discussion is
presented at ?buhring
) and the package employs
either hypergeo_gosper()
which uses iterative scheme due to
Gosper (Johansson et al. 2013), or the residue theorem if
The package comes with an extensive test suite in the tests/
directory. The tests fall into two main categories, firstly comparison
with either Maple or Mathematica output following Becken and Schmelcher (2000); and
secondly, verification of identities which appear in Abramowitz and Stegun (1965) as
elementary special cases. Consider, for example,
The left and right hand sides are given by eqn15.1.15a()
and eqn15.1.15b()
respectively which agree to numerical precision in
the test suite; but care must be taken with regard to the placing of
branch cuts. Further validation is provided by checking against known
analytical results. For example, it is known that
so, for example,
> hypergeo(2, 1, 2, -1/2)
[1] 0.66666666666667+0i
The hypergeo package offers direct numerical functionality to the R
user on the command line. The package is designed for use with R and
Figure 3 shows the package being used to
visualize
A second example is given from the author’s current work in game theory.
Consider a game in which a player is given
The two terms correspond to the final counter being removed from box
with R idiom:
> expected <- function(a, b, p) {
+ Re(
+ choose(a+b, b) * p^a * (1-p)^b *
+ (p * b/(1+a) * hypergeo(a+b+1, 2, a+2, p) +
+ (1-p) * a/(1+b) * hypergeo(a+b+1, 2, b+2, 1-p)))
+ }
Thus if
> c(expected(8, 2, 0.8), expected(9, 1, 0.8))
[1] 3.019899 1.921089
showing that the latter allocation is preferable in expectation.
Evaluation of the hypergeometric function is hard, as evidenced by the
extensive literature concerning its numerical
evaluation (Buhring 1987; Forrey 1997; Becken and Schmelcher 2000; Michel and Stoitsov 2008). The
hypergeo package is presented as a modular, R-centric implementation
with multiple evaluation methodologies, providing reasonably accurate
results over the complex plane and covering moderate real values of the
auxiliary parameters
NumericalMathematics, Optimization
This article is converted from a Legacy LaTeX article using the texor package. The pdf version is the official version. To report a problem with the html, refer to CONTRIBUTE on the R Journal homepage.
Hypergeometric2F1.pdf
, downloaded from http:/functions.wolfram.com/HypergeometricFunctions/.
Text and figures are licensed under Creative Commons Attribution CC BY 4.0. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".
For attribution, please cite this work as
Hankin, "Numerical Evaluation of the Gauss Hypergeometric Function with the hypergeo Package", The R Journal, 2015
BibTeX citation
@article{RJ-2015-022, author = {Hankin, Robin K. S.}, title = {Numerical Evaluation of the Gauss Hypergeometric Function with the hypergeo Package}, journal = {The R Journal}, year = {2015}, note = {https://rjournal.github.io/}, volume = {7}, issue = {2}, issn = {2073-4859}, pages = {81-88} }