I tried it. That is, in qhyper.c I replaced
p *= 1 - 64*DBL_EPSILON;
with
p *= 1 - 1000*DBL_EPSILON;
and recompiled. 2.3.0 alpha now passses make check-all and returns
sensible results on both of my FreeBSD 6.1 machines.
Rhyper <- 1:20
Phyper <- phyper (Rhyper, m = 40, n = 30, k = 20)
qhyper(Phyper, m = 40, n = 30, k = 20)
[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
Thanks very much for your guidance, this has been educational. Is
there anything else that I should check?
On Tue, Mar 28, 2006 at 01:53:09PM +0100, Prof Brian Ripley wrote:
qgamma.c comtains the line
p *= 1 - 64*DBL_EPSILON;
which seems to be too small nowadays. I find for example
p <- phyper(7, m = 40, n = 30, k = 20)
qhyper(p*(1+5e-16), m = 40, n = 30, k = 20)
qhyper(p*(1+6e-16), m = 40, n = 30, k = 20)
[1] 8
so it seems far more sensitive than the tolerance indicates is intended.
Now, phyper has been changed since that tolerance was added, and that
may be part of the explanation.
Can you try a larger tolerance, e.g 1000?
A better fix is probably not to use a tolerance but rather to do a
downwards search using phyper.
On Tue, 28 Mar 2006, Andrew Robinson wrote:
Phyper <- phyper(7, m = 40, n = 30, k = 20)
print(Phyper, digits=16)
[1] 0.01796062766370490
This is from 2.3.0:
Phyper <- phyper(7, m = 40, n = 30, k = 20)
print(Phyper, digits=16)
qhyper(Phyper, m = 40, n = 30, k = 20)
[1] 8
Then if I save Phyper from 2.2.1, and load it in 2.3.0,
qhyper(Phyper, m = 40, n = 30, k = 20)
[1] 7
so it appears that the difference between 0.01796062766370490 and
0.01796062766370491 is important here.
(2.3.0)
qhyper(0.01796062766370490, m = 40, n = 30, k = 20)
qhyper(0.01796062766370491, m = 40, n = 30, k = 20)
[1] 8
Comparing the full sample:
This is from 2.2.1:
1: 16 11 11 15 11 13 13 12 13 10 10 7 11 14 13 9 14 13 13 11
21:
Read 20 items
Phyper <- phyper (Rhyper, m = 40, n = 30, k = 20)
print(Phyper, digits=16)
[1] 0.99744478362220412 0.51295659016196715 0.51295659016196715
[4] 0.98680472393309626 0.51295659016196715 0.86627412777488222
[7] 0.86627412777488222 0.71494883441868140 0.86627412777488222
[10] 0.30864259597126753 0.30864259597126753 0.01796062766370490
[13] 0.51295659016196715 0.95139460528774533 0.86627412777488222
[16] 0.15132082044442866 0.95139460528774533 0.86627412777488222
[19] 0.86627412777488222 0.51295659016196715
This is from 2.3.0:
[1] 0.99744478362220401 0.51295659016196726 0.51295659016196726
[4] 0.98680472393309615 0.51295659016196726 0.86627412777488211
[7] 0.86627412777488211 0.71494883441868129 0.86627412777488211
[10] 0.30864259597126747 0.30864259597126747 0.01796062766370491
[13] 0.51295659016196726 0.95139460528774533 0.86627412777488211
[16] 0.15132082044442868 0.95139460528774533 0.86627412777488211
[19] 0.86627412777488211 0.51295659016196726
The 14's (elements 14 and 17) are identical, everything else is
slightly different.
I also got this interesting result from 2.2.1:
Rhyper <- 1:20
Phyper <- phyper (Rhyper, m = 40, n = 30, k = 20)
qhyper(Phyper, m = 40, n = 30, k = 20)
[1] 1 2 3 4 6 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
Repeating the operation for 2.3.0:
Rhyper <- 1:20
Phyper <- phyper (Rhyper, m = 40, n = 30, k = 20)
qhyper(Phyper, m = 40, n = 30, k = 20)
[1] 1 2 3 4 6 6 8 8 9 10 11 12 13 14 15 16 17 18 19 20
(I also tried on WinXP R 2.2.1, and got the expected 1:20 back.)
On Tue, Mar 28, 2006 at 12:02:11PM +0100, Prof Brian Ripley wrote:
Then we need you to dig in to find out why. I'd start by seeing if Phyper
had changed (either by printing it to 16dp or by saving it from 2.2.1 and
reloading into 2.3.0).
On Tue, 28 Mar 2006, Andrew Robinson wrote:
1: 16 11 11 15 11 13 13 12 13 10 10 7 11 14 13 9 14 13 13 11
21:
Read 20 items
Phyper <- phyper (Rhyper, m = 40, n = 30, k = 20)
qhyper(Phyper, m = 40, n = 30, k = 20)
[1] 16 11 11 15 11 13 13 12 13 10 10 8 11 14 13 9 14 13 13 11
The 12th element (8) differs from the input (7).
Phyper <- phyper (7, m = 40, n = 30, k = 20)
qhyper(Phyper, m = 40, n = 30, k = 20)
[1] 8
If I do this using 2.2.1 then the input and the output are identical.
On Tue, Mar 28, 2006 at 10:27:08AM +0100, Prof Brian Ripley wrote:
On Tue, 28 Mar 2006, Andrew Robinson wrote:
You're welcome. You are correct. d-p-q-r-tests.Rout.fail
shows:
All.eq(Rhyper, qhyper (Phyper, m = 40, n = 30, k = 20))
[1] "Mean scaled difference: 0.08333333"
Yes, please run the lines below, e.g.
Rhyper <- scan()
16 11 11 15 11 13 13 12 13 10 10 7 11 14 13 9 14 13 13 11
Phyper <- phyper (Rhyper, m = 40, n = 30, k = 20)
qhyper(Phyper, m = 40, n = 30, k = 20)
and tell us what answer you get.
Let me know if/how I can further assist.
Andrew
On Tue, Mar 28, 2006 at 09:03:48AM +0100, Prof Brian Ripley wrote:
Thanks for checking.
Please look in d-p-q-r-tests.Rout.fail and see what immediately
preceeds
the line
[1] "Mean scaled difference: 0.08333333"
Some experimentation suggests it is
All.eq(Rhyper, qhyper (Phyper, m = 40, n = 30, k = 20))
If so, we have
Rhyper <- scan()
16 11 11 15 11 13 13 12 13 10 10 7 11 14 13 9 14 13 13 11
Phyper <- phyper (Rhyper, m = 40, n = 30, k = 20)
and those have been checked. So the error would appear to be in
qhyper(Phyper, m = 40, n = 30, k = 20)
and indeed a mean scaled difference of 1/12 is plausible, since the
mean
of Rhyper is 12. So I deduce that your platform has a problem in
qhyper,
but please cross-check.
If so, this is strange as the only recent change to qhyper.c (or
things
I
can see it uses such as lfastchoose) is cosmetic.
Can you confirm the diagnosis is correct so far?
On Tue, 28 Mar 2006, Andrew Robinson wrote:
Hi Developers,
The alpha, compiles successfully, but it is failing make check-all
(on
two seperate machines, both FreeBSD 6.1).
Here is the version string:
platform i386-unknown-freebsd6.1
arch i386
os freebsd6.1
system i386, freebsd6.1
status alpha
major 2
minor 3.0
year 2006
month 03
day 27
svn rev 37584
language R
version.string Version 2.3.0 alpha (2006-03-27 r37584)
Here is the error message from make check-all
comparing 'd-p-q-r-tests.Rout' to './d-p-q-r-tests.Rout.save'
...706c706
< [1] "Mean scaled difference: 0.08333333"
---
gmake[3]: *** [d-p-q-r-tests.Rout] Error 1
gmake[3]: Leaving directory `/usr/local/beta/R-alpha/tests'
gmake[2]: *** [test-Specific] Error 2
gmake[2]: Leaving directory `/usr/local/beta/R-alpha/tests'
gmake[1]: *** [test-all-basics] Error 1
gmake[1]: Leaving directory `/usr/local/beta/R-alpha/tests'
gmake: *** [check-all] Error 2
I have checked d-p-q-r-tests.Rout.fail for any obvious problems - I
found some warnings, viz.
pgamma(1,Inf,scale=Inf) == 0
## Also pgamma(Inf,Inf) == 1 for which NaN was slightly more
all(is.nan(c(pgamma(Inf, 1,scale=Inf),
+ pgamma(Inf,Inf,scale=Inf))))
[1] TRUE
Warning messages:
1: NaNs produced in: pgamma(q, shape, scale, lower.tail, log.p)
2: NaNs produced in: pgamma(q, shape, scale, lower.tail, log.p)
x0 <- -2 * 10^-c(22,10,7,5)
stopifnot(pbinom(x0, size = 3, prob = 0.1) == 0,
+ dbinom(x0, 3, 0.1) == 0) # d*() warns about non-integer
Warning messages:
1: non-integer x = -0.000000
2: non-integer x = -0.000020
## very small negatives were rounded to 0 in R 2.2.1 and earlier
I hope that this is helpful. Thanks are due to Peter Dalgaard for
guidance. So, thanks Peter :).
Cheers
Andrew
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595