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]]
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
