Certainly. I assume that you mean qhyper.c rather than qgamma.c. Can you please be more explicit as to how I can try your suggestion. Do you mean that I should try:
p *= 1 - 1000*DBL_EPSILON; ? Thanks. 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) > [1] 7 > >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: > > >This is from 2.2.1: > > > >>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) > >[1] 0.01796062766370491 > > > >>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) > >[1] 7 > >>qhyper(0.01796062766370491, m = 40, n = 30, k = 20) > >[1] 8 > > > > > > > >Comparing the full sample: > > > >This is from 2.2.1: > > > >>Rhyper <- scan() > >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: > > > >>print(Phyper, digits=16) > >[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: > >> > >>>I get: > >>> > >>>>Rhyper <- scan() > >>>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" > >>>>>>>--- > >>>>>>>>[1] TRUE > >>>>>>>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 > >>>>>>>[1] TRUE > >>>>>>>>## Also pgamma(Inf,Inf) == 1 for which NaN was slightly more > >>>>>>>appropriate > >>>>>>>>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, [EMAIL PROTECTED] > >>>>>>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, [EMAIL PROTECTED] > >>>>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, [EMAIL PROTECTED] > >>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, [EMAIL PROTECTED] > 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 -- Andrew Robinson Department of Mathematics and Statistics Tel: +61-3-8344-9763 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 Email: [EMAIL PROTECTED] http://www.ms.unimelb.edu.au ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
