[Rd] pt inaccurate when x is close to 0 (PR#9945)

[jerry.lewis at biogenidec.com]

While I agree that the reported results from Mathematica have only 10-13 
correct digits, that does not mean that pt() in R is any better for these 
calculations.  For instance the following three calculations are 
mathematically equivalent, but pt() disagrees at the 13th figure in R 
v2.6.2
  pt(1e-4,13)
  pf(1e-8,1,13)/2+0.5
  pbeta(1/(1+13/1e-8),.5,6.5)/2+0.5
Using 1/(1+n/x^2) and reversing the shape parameters avoids the accuracy 
loss that Thomas Lumley was concerned about.

The following reference values for pt(1e-4,1:100) were computed in Maple 
to 400 figures and then rounded to 17.  To 17 figures, there was no 
difference in the results whether I used the exact value of 10^(-4) or the 
exact binary double precision approximation of 7378697629483821*2^(-66). 
Compared to these reference values, log relative errors (LRE, essentially 
the number of correct figures) for pt ranged from 10.7 to 13.4, while
  pf(1e-8,1,df)/2+0.5
and
  pbeta(1/(1+df/1e-8),.5,df/2)/2+0.5
agreed exactly with Maple, except for 21 and 90 df, where the LRE was 
15.65.  Clearly it would be quite easy to improve pt() for small values of 
x.


Jerry W. Lewis
Lead Statistician
Biogen Idec
Cambridge, MA  USA


z<-c(0.50003183098851228, 0.50003535533897094, 0.50003675525961311, 
0.50003749999992188, 0.50003796066890633, 0.50003827327715657, 
0.50003849914500990, 0.50003866990202363, 0.50003880349081531, 
0.50003891083832527, 0.50003899897563476, 0.50003907263045176, 
0.50003913509812855, 0.50003918874550927, 0.50003923531621426, 
0.50003927612297732, 0.50003931217293135, 0.50003934425160310, 
0.50003937298066018, 0.50003939885850220, 0.50003942228935944, 
0.50003944360452316, 0.50003946307808180, 0.50003948093875281, 
0.50003949737889518, 0.50003951256145633, 0.50003952662538506, 
0.50003953968989187, 0.50003955185783298, 0.50003956321842127, 
0.50003957384941549, 0.50003958381890103, 0.50003959318674852, 
0.50003960200581637, 0.50003961032294799, 0.50003961817980360, 
0.50003962561355779, 0.50003963265748728, 0.50003963934146876, 
0.50003964569240221, 0.50003965173457262, 0.50003965748996017, 
0.50003966297850729, 0.50003966821834945, 0.50003967322601522, 
0.50003967801660046, 0.50003968260392016, 0.50003968700064152, 
0.50003969121840071, 0.50003969526790563, 0.50003969915902671, 
0.50003970290087718, 0.50003970650188431, 0.50003970996985278, 
0.50003971331202108, 0.50003971653511200, 0.50003971964537768, 
0.50003972264864014, 0.50003972555032763, 0.50003972835550734, 
0.50003973106891495, 0.50003973369498129, 0.50003973623785651, 
0.50003973870143192, 0.50003974108935987, 0.50003974340507184, 
0.50003974565179484, 0.50003974783256646, 0.50003974995024854, 
0.50003975200753972, 0.50003975400698688, 0.50003975595099569, 
0.50003975784184024, 0.50003975968167193, 0.50003976147252762, 
0.50003976321633717, 0.50003976491493036, 0.50003976657004330, 
0.50003976818332434, 0.50003976975633955, 0.50003977129057782, 
0.50003977278745549, 0.50003977424832079, 0.50003977567445783, 
0.50003977706709040, 0.50003977842738546, 0.50003977975645637, 
0.50003978105536601, 0.50003978232512953, 0.50003978356671702, 
0.50003978478105603, 0.50003978596903380, 0.50003978713149950, 
0.50003978826926618, 0.50003978938311274, 0.50003979047378564, 
0.50003979154200064, 0.50003979258844427, 0.50003979361377541, 
0.50003979461862660)
	[[alternative HTML version deleted]]