-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1
<<insert bug report here>>
Reproduced on Debian and Windows ...
On 2.4.x if you execute
set.seed(12345)
t.1 <- rt(n = 1000, df = 20)
set.seed(12345)
t.2 <- rt(n = 1000, df = 20, ncp = 0)
all.equal(t.1, t.2) ## Not close to true
This appears to be due to the fact that in 2.4.x rt is now
rt
function (n, df, ncp = 0)
{
if (missing(ncp))
.Internal(rt(n, df))
else rnorm(n, ncp)/sqrt(rchisq(n, df)/df)
}
<environment: namespace:stats>
Whereas in 2.3.1 rt() is verified to work as expected when someone
(redundantly) types ncp = 0
rt
function (n, df, ncp = 0)
{
if (ncp == 0)
.Internal(rt(n, df))
else rnorm(n, ncp)/(rchisq(n, df)/sqrt(df))
}
<environment: namespace:stats>
The only interface difference is missing(ncp) vs. ncp == 0 .
- --please do not edit the information below--
Version:
platform = i486-pc-linux-gnu
arch = i486
os = linux-gnu
system = i486, linux-gnu
status major = 2
minor = 4.1
year = 2006
month = 12
day = 18
svn rev = 40228
language = R
version.string = R version 2.4.1 (2006-12-18)
Locale:
LC_CTYPE=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8;LC_COLLATE=en_US.UTF-8;LC_MONETARY=en_US.UTF-8;LC_MESSAGES=en_US.UTF-8;LC_PAPER=en_US.UTF-8;LC_NAME=en_US.UTF-8;LC_ADDRESS=en_US.UTF-8;LC_TELEPHONE=en_US.UTF-8;LC_MEASUREMENT=en_US.UTF-8;LC_IDENTIFICATION=en_US.UTF-8
Search Path:
.GlobalEnv, package:car, package:rkward, package:stats,
package:graphics, package:grDevices, package:utils, package:datasets,
package:methods, Autoloads, package:base
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.6 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org
iD8DBQFFiF09zQDSXIcN85kRAhK4AKCBmPNkGszTzFdALfiGgW1LfNIOTwCfZ8RM
0gZjer8+rnmOPDFKZXBbFQw=ff85
-----END PGP SIGNATURE-----
>>>>> "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
Why do you consider this might be a bug in R rather than in your expectations? Your 'expected' is not what others have expected .... In 2.4.1 rt(ncp) operates continuously as ncp is varied through 0. That is different from saying that you want a central t (omitting ncp). BTW, reports to R-bugs should *never* include a question mark: see the R FAQ. A question on R-devel (or even R-help) would have been reasonable. On Tue, 19 Dec 2006, goodrich at fas.harvard.edu wrote:> -----BEGIN PGP SIGNED MESSAGE----- > Hash: SHA1 > > <<insert bug report here>> > Reproduced on Debian and Windows ... > > On 2.4.x if you execute > > set.seed(12345) > t.1 <- rt(n = 1000, df = 20) > > set.seed(12345) > t.2 <- rt(n = 1000, df = 20, ncp = 0) > > all.equal(t.1, t.2) ## Not close to true > > This appears to be due to the fact that in 2.4.x rt is now > > rt > function (n, df, ncp = 0) > { > if (missing(ncp)) > .Internal(rt(n, df)) > else rnorm(n, ncp)/sqrt(rchisq(n, df)/df) > } > <environment: namespace:stats> > > Whereas in 2.3.1 rt() is verified to work as expected when someone > (redundantly) types ncp = 0 > > rt > function (n, df, ncp = 0) > { > if (ncp == 0) > .Internal(rt(n, df)) > else rnorm(n, ncp)/(rchisq(n, df)/sqrt(df)) > } > <environment: namespace:stats> > > The only interface difference is missing(ncp) vs. ncp == 0 . > > - --please do not edit the information below-- > > Version: > platform = i486-pc-linux-gnu > arch = i486 > os = linux-gnu > system = i486, linux-gnu > status > major = 2 > minor = 4.1 > year = 2006 > month = 12 > day = 18 > svn rev = 40228 > language = R > version.string = R version 2.4.1 (2006-12-18) > > Locale: > LC_CTYPE=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8;LC_COLLATE=en_US.UTF-8;LC_MONETARY=en_US.UTF-8;LC_MESSAGES=en_US.UTF-8;LC_PAPER=en_US.UTF-8;LC_NAME=en_US.UTF-8;LC_ADDRESS=en_US.UTF-8;LC_TELEPHONE=en_US.UTF-8;LC_MEASUREMENT=en_US.UTF-8;LC_IDENTIFICATION=en_US.UTF-8 > > Search Path: > .GlobalEnv, package:car, package:rkward, package:stats, > package:graphics, package:grDevices, package:utils, package:datasets, > package:methods, Autoloads, package:base > -----BEGIN PGP SIGNATURE----- > Version: GnuPG v1.4.6 (GNU/Linux) > Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org > > iD8DBQFFiF09zQDSXIcN85kRAhK4AKCBmPNkGszTzFdALfiGgW1LfNIOTwCfZ8RM > 0gZjer8+rnmOPDFKZXBbFQw> =ff85 > -----END PGP SIGNATURE----- > > ______________________________________________ > 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
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 Prof Brian Ripley wrote:> Why do you consider this might be a bug in R rather than in your > expectations? Your 'expected' is not what others have expected ....Thank you for this response. I did consider at the time that this unexpected result could be deliberate, but help(rt) on 2.4.1 seemed to confirm my expectations. Martin Maechler noted on the first response to this bug report that help(rt) was not updated to describe the new functionality of rt().> In 2.4.1 rt(ncp) operates continuously as ncp is varied through 0. > That is different from saying that you want a central t (omitting ncp).I was not expecting random draws from a central t distribution to be the same as random draws from a limiting case of a non-central t distribution. I was expecting that redundantly typing the default arguments of a function would produce the same result as omitting the default arguments. Are there other functions with this discrepancy? In other words, I do not understand why rt() continues to default to ncp = 0, as opposed to no default argument for ncp. It would seem that args(rt) currently implies that the draws are to be taken from a limiting case of non-central t distribution, when this is not accurate if the ncp argument is omitted. Thanks, Ben -----BEGIN PGP SIGNATURE----- Version: GnuPG v1.4.6 (GNU/Linux) Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org iD8DBQFFid58zQDSXIcN85kRAhRVAJ0WNf2fMVcPZk4rfxw8qXO272Bs8ACfdVAt a5OaqyObzqkjL8a5sCooj6w=/fPt -----END PGP SIGNATURE-----