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