From maechler at stat.math.ethz.ch  Wed Dec 20 08:50:49 2006
From: maechler at stat.math.ethz.ch (maechler at stat.math.ethz.ch)
Date: Wed, 20 Dec 2006 08:50:49 +0100 (CET)
Subject: [Rd] Bug in rt() ? (PR#9422)
Message-ID: <20061220075049.856F74C9F8@slim.kubism.ku.dk>

>>>>> "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



[Rd] Bug in rt() ? (PR#9422)
maechler at stat.math.ethz.ch