From ripley at stats.ox.ac.uk  Fri Feb 15 17:06:23 2008
From: ripley at stats.ox.ac.uk (Prof Brian Ripley)
Date: Fri, 15 Feb 2008 16:06:23 +0000 (GMT)
Subject: [Rd] pt inaccurate when x is close to 0 (PR#9945)
In-Reply-To: <20080215152008.0C4D8283462C@mail.pubhealth.ku.dk>
References: <20080215152008.0C4D8283462C@mail.pubhealth.ku.dk>
Message-ID: <alpine.LFD.1.00.0802151603390.4041@gannet.stats.ox.ac.uk>

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



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