Re: Numerical accuracy of Hypergeometric2F1

• To: mathgroup at smc.vnet.net
• Subject: [mg55754] Re: Numerical accuracy of Hypergeometric2F1
• From: Maxim <ab_def at prontomail.com>
• Date: Tue, 5 Apr 2005 03:22:02 -0400 (EDT)
• References: <d2qj7t\$p0o\$1@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

```On Mon, 4 Apr 2005 05:27:25 +0000 (UTC), Christos Argyropoulos M.D.
<chrisarg at fuse.net> wrote:

>
> Hi,
> Recently I ran into a problem in applied statistics which required
> evaluation of a specific Hypergeometric2F1 functional i.e.
> Hypergeometric2F1[k+1/2,1,3/2,x] where Element[k,Integers], k>0 ,
> Element[x, Reals] and  0<=x<=1.
> It appears that for "large" values of k , Mathematica returns the wrong
> If one re-writes the hypergeometric as a polynomial (using the distant
> neighbor relation  http://functions.wolfram.com/07.23.17.0007.01), get
> rid of the ratio of Pochhamer functions and then examine the numerical
> values given by the two formulations diverge.
> \!\(fpol[x_, k_] :=
>     Module[{}, coeff = Table[1, {k + 1}];
>       For[i = 2, i \[LessEqual] k + 1, \(i++\),
>         coeff[\([i]\)] =
>           coeff[\([i -
>                   1]\)]\ \(\(-k\) + i - 1\)\/\(i - k - 1/2\)]; \
> {Hypergeometric2F1[k + 1/2, 1, 3/2,
>           x], \(1\/\(2  k -
>                 1\)\) \(\[Sum]\+\(i = 1\)\%\(k + 1\)\((\((1 -
> x)\)\^\(-i\)\ \
> coeff[\([i]\)])\)\)}]\)
>
> In[4]:=
> SetPrecision[fpol[0.7,100],16]
> Out[4]=
> \!\({\(-2.5109638442451095387366941924357`15.9546*^56\),
>     2.0578481017803705921545468593643`15.9546*^51}\)
>
> but
>
> In[5]:=
> SetPrecision[fpol[0.7,10],16]
> Out[5]=
> {57439.31727367684,57439.31727367681}
>
> or graphically ...
> ListPlot3D[Table[Apply[(#1/#2)&,fpol[x,k]],{x,0.,1.,.1},{k,1,100}]]
>
> The problems do not seem to be caused by numerical errors secondary to
> direct evaluation of high - order polynomials in the equivalent
> expression since, using the Horner formulation for the polynomial gives
> essentially identical results
>
> \!\(fhorn[x_, k_] :=
>     Module[{}, coeff = Table[1, {k + 1}];
>       For[i = 2, i \[LessEqual] k + 1, \(i++\),
>         coeff[\([i]\)] =
>           coeff[\([i -
>                   1]\)]\ \(\(-k\) + i - 1\)\/\(i - k - 1/2\)]; \
> {Hypergeometric2F1[k + 1/2, 1, 3/2, x], \(1\/\(2  k - 1\)\)
>           Module[{y = 1/\((1 - x)\), res = 1}, res = 0;
>             For[i = k + 1, i \[GreaterEqual] \ 1\ , \(i--\),
>               res = \((res + coeff[\([i]\)])\)\ y\ ]; res]}]\)
>
> In[6]:=
> SetPrecision[fhorn[0.7,10],16]
> Out[6]=
> {57439.31727367684,57439.31727367682}
>
> and
>
> In[7]:=
> SetPrecision[fhorn[0.7,100],16]
> Out[7]=
> \!\({\(-2.5109638442451095387366941924357`15.9546*^56\),
>     2.0578481017803822228995099773782`15.9546*^51}\)
>
> and graphically :
> ListPlot3D[Table[Apply[(#1/#2)&,fhorn[x,k]],{x,0.,1.,.1},{k,1,100}]]
>
> Has anyone else run into similar problems with the Hypergeometric2F1 for
> large values of one of the parameters? Is there a quick and dirty way of
> knowing when to mistruct the hypergeometric numerical answer returned by
> Mathematica ?
>
> Christos Argyropoulos MD PhD
> University of Cincinnati,
> College of Medicine
>

First note that SetPrecision isn't very useful here. You compute
Hypergeometric2F1 with machine precision (Mathematica sees the machine
number 0.7 and then all the subsequent computations are performed using
hardware floating-point numbers) and only after that tell SetPrecision to
convert the result to an arbitrary-precision number, padding its binary
representation with zeros if necessary.

Machine arithmetic provides no precision control. At best Mathematica can
look at the arguments and decide which numerical method is more suitable
for the given domain of their values; perhaps this decision process can be
improved, but in general machine-precision computations certainly can give
wrong results.

The problems begin when we turn to significance arithmetic (the outputs
are for version 5.1.0):

In[1]:=
Hypergeometric2F1[k + 1/2, 1, 3/2, 7/10] /. k -> 1000`20

Out[1]=
-3.08039205255934313614253847`16.439170942121546*^720

In[2]:=
Hypergeometric2F1[k + 1/2, 1, 3/2, 7/10] /. k -> 1000 // N[#, 20]&

Out[2]=
2.53393510249899035708859330734`20.*^521

It is easy to verify the result because if we evaluate Hypergeometric2F1
with exact values for k and x (In[2]), then we obtain a rational number,
numerical evaluation of which doesn't present a problem. Thus we can see
that no digits in Out[1] are correct, even though Mathematica 'guarantees'
us 16 significant digits.

So in this sense your observation is correct: Mathematica's significance
arithmetic doesn't really help here, especially when the parameters a,b,c
are large or close to negative integer values. Another case that presents
difficulties for numerical evaluation is the computation of the
derivatives with respect to the parameters. Then we just have to resort to
the good old method of reevaluating the input at higher precision and
comparing the results.

Maxim Rytin
m.r at inbox.ru

```

• Prev by Date: Having trouble with substitution tile at higher iteration levels--> takes forever!
• Next by Date: Recommendations for a programming book?
• Previous by thread: Re: Numerical accuracy of Hypergeometric2F1
• Next by thread: Re: Re: Numerical accuracy of Hypergeometric2F1