R devel - Jan 2008 - chi-squared with zero df (PR#10551)

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


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

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

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

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

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

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

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

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

Not on my Linux computers; e.g.,

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

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

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

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

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

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

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

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

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

i.e. "numerically 0"

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

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

I'll investigate.

Martin Maechler, ETH Zurich

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