Re: [Rd] Bug in rt() ? (PR#9422)

From: <maechler_at_stat.math.ethz.ch>
Date: Wed 20 Dec 2006 - 07:50:49 GMT


>>>>> "goodrich" == goodrich <goodrich@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



R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Thu Dec 21 19:03:57 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Thu 21 Dec 2006 - 08:31:03 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.