jerry.lewis at biogenidec.com
2008-Feb-15 15:20 UTC
[Rd] pt inaccurate when x is close to 0 (PR#9945)
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]]
Please do as the FAQ asked you to, and check the NEWS file at https://svn.R-project.org/R/trunk/NEWS R 2.6.2 is not the '_next_ release of R', and that pre-release behaves differently as documented in the NEWS file. On Fri, 15 Feb 2008, jerry.lewis at biogenidec.com wrote:> 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 at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel >-- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595