[Rd] Noncentral dt() with tiny 'x' values (PR#8874)

[maechler at stat.math.ethz.ch]

>>>>> "MikeMer" == meredith  <meredith at easynet.co.uk>
>>>>>     on Thu, 18 May 2006 03:52:51 +0200 (CEST) writes:

    MikeMer> Full_Name: Mike Meredith
    MikeMer> Version: 2.3.0
    MikeMer> OS: WinXP SP2
    MikeMer> Submission from: (NULL) (210.195.228.29)



    MikeMer> Using dt() with a non-centrality parameter and near-zero values for 'x' results
    MikeMer> in erratic output. Try this:

    MikeMer> tst <- c(1e-12, 1e-13, 1e-14, 1e-15, 1e-16, 1e-17, 0)
    MikeMer> dt(tst,16,1)

    MikeMer> I get:  0.2381019 0.2385462 0.2296557 0.1851817 0.6288373 3.8163916 (!!)
    MikeMer> 0.2382217

I get quite different values (several '0', BTW), which just
confirms the erratic nature.

As often, plots give even a clearer picture:

if(!require(sfsmisc)) # pkg 'sfsmisc'
    lseq <- function (from, to, length) 
         2^seq(log2(from), log2(to), length.out = length)
x <- lseq(1e-3, 1e-33, length= 301)

plot(x, dt(x, df=16, ncp=1),   type = "o", cex=.5, log = "x")
plot(x, dt(x, df=16, ncp=0.1), type = "o", cex=.5, log = "x")
plot(x, dt(x, df= 3, ncp=0.1), type = "o", cex=.5, log = "x")



    MikeMer> The 0.238 values are okay, the others nonsense, and
    MikeMer> they cause confusing spikes on plots of dt() vs 'x'
    MikeMer> if 'x' happens to include tiny values. (Other
    MikeMer> values of df and ncp also malfunction, but not all
    MikeMer> give results out by an order of magnitude!)

I think almost all do, once you start looking at plots like the above.

    MikeMer> I'm using the work-around dt(round(x,10),...), but
    MikeMer> dt() should really take care of this itself.

or actually rather do something more smart; the cutoff at 1e-10
is quite crude.

Note that this is not a new bug at all; but rather as old as
we have dt(*, ncp= .) in R.

Thanks for reporting it!
Martin Maechler, ETH Zurich