[Date Index] [Thread Index] [Author Index]

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
> answer .
> 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