[Rd] Bug in rt() ? (PR#9422)
maechler at stat.math.ethz.ch
maechler at stat.math.ethz.ch
Wed Dec 20 08:50:49 CET 2006
>>>>> "goodrich" == goodrich <goodrich at fas.harvard.edu>
>>>>> on Tue, 19 Dec 2006 22:45:33 +0100 (CET) writes:
goodrich> Reproduced on Debian and Windows ...
goodrich> On 2.4.x if you execute
goodrich> set.seed(12345)
goodrich> t.1 <- rt(n = 1000, df = 20)
goodrich> set.seed(12345)
goodrich> t.2 <- rt(n = 1000, df = 20, ncp = 0)
goodrich> all.equal(t.1, t.2) ## Not close to true
goodrich> This appears to be due to the fact that in 2.4.x rt is now
goodrich> rt
goodrich> function (n, df, ncp = 0)
goodrich> {
goodrich> if (missing(ncp))
goodrich> .Internal(rt(n, df))
goodrich> else rnorm(n, ncp)/sqrt(rchisq(n, df)/df)
goodrich> }
goodrich> <environment: namespace:stats>
goodrich> Whereas in 2.3.1 rt() is verified to work as expected when someone
goodrich> (redundantly) types ncp = 0
goodrich> rt
goodrich> function (n, df, ncp = 0)
goodrich> {
goodrich> if (ncp == 0)
goodrich> .Internal(rt(n, df))
goodrich> else rnorm(n, ncp)/(rchisq(n, df)/sqrt(df))
goodrich> }
goodrich> <environment: namespace:stats>
goodrich> The only interface difference is missing(ncp) vs. ncp == 0 .
Your analysis is correct.
But the question mark in your subject line is too: This is no
bug, but quite on purpose: rt() has been changed to be in line with
qt() and pt().
The idea is that if you specify 'ncp', we give you behavior which
is continuous (in the mathematical sense) in 'ncp' for ncp --> 0.
Hence, when ncp is specified {aka !missing(.)}, the algorithms
for ncp >=0 are applied; if it's not specified, the algorithms
of the central t-distributions are used.
However, your report does point to a ``bug'' :
We should document the above on the Tdist help page.
So, after all, thank you for the report!
Martin Maechler, ETH Zurich
More information about the R-devel
mailing list