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

maechler at stat.math.ethz.ch maechler at stat.math.ethz.ch
Mon Jan 7 09:50:15 CET 2008


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



More information about the R-devel mailing list