[Rd] chi-squared with zero df (PR#10551)

[maechler at stat.math.ethz.ch]

>>>>> "MM" == Martin Maechler <maechler at stat.math.ethz.ch>
>>>>>     on Mon,  7 Jan 2008 09:50:15 +0100 (CET) writes:

>>>>> "JL" == Jerry Lewis <Jerry.Lewis at biogenidec.com>
>>>>>     on Mon,  7 Jan 2008 05:20:23 +0100 (CET) writes:

    JL> Full_Name: Jerry W. Lewis
    JL> Version: 2.6.1
    JL> OS: Windows XP Professional
    JL> Submission from: (NULL) (24.147.191.250)

    JL> pchisq(0,0,ncp=lambda) returns 0 instead of exp(-lambda/2)

    MM> Yes. As we know,  chi-square(df=0, ncp=lambda)   has a point mass of
    MM> exp(-lambda/2) at 0.

    MM> Hence pchisq() has a corresponding jump there; strictly, it's a
    MM> matter of definition (left- / right- continuity, etc) what
    MM> pchisq() should be there, but indeed, the usual definition --
    MM> which we also follow for the discrete distributions -- 
    MM> is "cadlag" (continue à droite, limite à gauche), 
    MM> and you are right.

    MM> That's easily fixed.

    JL> pchisq(x,0,ncp=lambda) returns NaN 

    JL> instead of exp(-lambda/2)*(1 + SUM_{r=0}^infty
    JL> ((lambda/2)^r / r!) pchisq(x, df + 2r))

    MM> Not on my Linux computers; e.g.,

    >> pchisq(0:10, 0,1)
    MM> [1] 0.0000000 0.7328798 0.8193100 0.8781745 0.9181077 0.9451020 0.9632911
    MM> [8] 0.9755110 0.9836985 0.9891706 0.9928194

    MM> but I can see the problem on Windows (Server 2003 R2, x64 edition),
    MM> where I get 

    >> pchisq(0:10, 0,1)
    MM> [1] 0 NaN NaN ... NaN

    MM> aah, and I also see it on a 32-bit Linux computer.

and actually did see it on all computers; I was wrong above.

It's been a consequence of an "improvement" introduced for R 2.3.0;
unfortunately that made the (df=0, ncp < 80) wrongly return NaN.

Now fixed in both R-patched and R-devel.

Thank you once more for the bug report!
Martin


    JL> qchisq(.7,0,ncp=1) returns 1.712252 instead of 0.701297103

    MM> On my 64-bit Linux computer I get the correct result where as
    MM> the above Windows and our 32-bit Linux give *different* (but
    MM> wrong) results.

    JL> qchisq(exp(-1/2),0,ncp=1) returns 1.238938 instead of 0

    MM> Not on my 64-bit Linux where e.g.,

    >> lam <- 1:10; qchisq(exp(-lam), 0, ncp=lam)
    MM> [1] 2.225074e-308 2.225074e-308 2.225074e-308 2.225074e-308 2.225074e-308
    MM> [6] 2.225074e-308 2.225074e-308 2.225074e-308 2.225074e-308 2.225074e-308

    MM> i.e. "numerically 0"

    MM> {Note that I've known about problems with our non-central chisq
    MM> algorithms, but these were all of very extreme cases...}

    MM> In summary, we seem to have an issue with our algorithms failing
    MM> on some platforms.

    MM> I'll investigate.

    MM> Martin Maechler, ETH Zurich