>>>>> "JL" == Jerry Lewis <[EMAIL PROTECTED]>
>>>>> 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
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
