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

[maechler at stat.math.ethz.ch]

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