Christian.Stratowa@vie.boehringer-ingelheim.com
2003-Feb-14 11:29 UTC
[R] FW: [Fwd: Re: [S] Exact p-values]
Dear all Just for fun, I have just downloaded the paper mentioned below and checked it with R-1.6.1. Everything is ok with exception of Table 2b, where I get always 1 instead of 0.5:> pbinom(1e15,2e15,0.5)[1] 1 Which value should be correct? Best regards Christian Stratowa =============================================Christian Stratowa, PhD Boehringer Ingelheim Austria Dept NCE Lead Discovery - Bioinformatics Dr. Boehringergasse 5-11 A-1121 Vienna, Austria Tel.: ++43-1-80105-2470 Fax: ++43-1-80105-2683 email: christian.stratowa at vie.boehringer-ingelheim.com> -------- Original Message -------- > Subject: Re: [S] Exact p-values > Date: Thu, 13 Feb 2003 18:31:38 +0100 > From: "Rau, Roland" <Rau at demogr.mpg.de> > To: 'Spencer Graves' <spencer.graves at PDF.COM>, Jose Mar?a Fedriani > Laffitte <fedriani at ebd.csic.es> > CC: s-news at lists.biostat.wustl.edu > > Dear all, > > in relation to your question, the following working paper of Leo Knuesel, > University of Munich, might be of interest: > "On the Accuracy of Statistical Distributions in S-Plus for Windows > (1999)" > You can download the paper from (pdf-Format, 45k): > http://www.stat.uni-muenchen.de/~knuesel/elv/accuracy.html > > Best, > Roland > > > -----Original Message----- > > From: Spencer Graves [SMTP:spencer.graves at PDF.COM] > > Sent: Thursday, February 13, 2003 6:12 PM > > To: Jose Mar?a Fedriani Laffitte > > Cc: s-news at lists.biostat.wustl.edu > > Subject: Re: [S] Exact p-values > > > > > > Try ( 1-pchisq(29.8, df=1)): With S-Plus 6.1, I got 4.78992e-008. > > > > By the way, the distribtion functions in R have more > arguments. > > For example, pchisq(29.8, df=1, lower.tail=F) produces the same > > answer, and pchisq(29.8, df=1, lower.tail=F, log=T) produces its > natural > > logarithm. Also, pchisq, dchisq, qchisq, and rchisq in R all have an > > "ncp" noncentrality parameter argument; only pchisq has such in S-Plus > > 6.1. Similarly, none of the Student's t functions in S-Plus have a > > non-centralitity parameter; in R, pt has an argument ncp, and from > this > > one can easily program ncp for dt, qt and rt. Also, the distribution > > functions in the current release of S-Plus are known to have problems. > > For example, pt(-1, Inf) = 0.5 in S-Plus 6.1, but 0.159 in R; > clearly, > > S-Plus gives a wrong answer without warning. > > > > Best Wishes, > > Spencer Graves > > > > Jose Mar?a Fedriani Laffitte wrote: > > > > >Dear all, > > > > > > I want to get the exact p-values, on 1 degree of freedom, for an > > array > > >of chi-square values. When my chi-square values are equal or lower > than > > >29.7, I get the exact associated p-values. Thus, for instance: > > > > > > > > > > > >>pchisq(29.7, df=1) > > >> > > >> > > >[1] 0.9999999 > > > > > >However, when my chi-square values are greater or equal to 29.8 what I > > get > > >is: > > > > > > > > > > > >>pchisq(29.8, df=1) > > >> > > >> > > >[1] 1 > > > > > > > > > Could anyone tell me how to fix this trivial issue? Very > grateful, > > Jose > > >M. Fedriani > > > > > >**************************************** > > >Jose M? Fedriani Laffitte > > >Estacion Biologica de Donana (CSIC) > > >Avda. M? Luisa s/n > > >41013-Sevilla > > >Spain > > >Tel. +34-954232340 > > >Fax +34-954621125 > > >http://ebd.csic.es > > > > > >-------------------------------------------------------------------- > > >This message was distributed by s-news at lists.biostat.wustl.edu. To > > >...(s-news-request clipped)...> > > > > > > > > > > > -------------------------------------------------------------------- > > This message was distributed by s-news at lists.biostat.wustl.edu. To > > ...(s-news-request clipped)...> -------------------------------------------------------------------- > This message was distributed by s-news at lists.biostat.wustl.edu. To > ...(s-news-request clipped)...> >
To understand the correct answer, you need to understand the following: > pbinom(1, 2, .5) [1] 0.75 This is the binomial cumulative distribution function. *** pbinom(0, 2, .5) = 0.25 *** pbinom(1, 2, .5) = 0.75 = 0.25 + 0.5 *** pbinom(2, 2, .5) = 1 However, pbinom(1e15, 2e15, .5) is a computational challenge. Standard numerical algorithms often fail in situations like this. The code should test for such cases and use more numerically stable "approximations" in place of the "exact" algorithms. The standard deviation for a binomial is sqrt(p*(1-p)/n) = 0.5/sqrt(2e15), which is roughly 1e-8 in your case. I get the following from both S-Plus and R: > pbinom(1e5+c(-1, 0, 1), 2e5, .5) [1] 0.4991079 0.5008921 0.5026762 For the problem you cite, the correct answer should be 0.5 to about 8 significant digits. Instead, I get 1 from R (as you did) and the following from S-Plus: > pbinom(1e15,2e15,0.5) [1] 0.7411209 Both give wrong answers without warning, though in this case, S-Plus is closer. Answer the question? Spencer Graves ##################### Christian.Stratowa at vie.boehringer-ingelheim.com wrote:> Dear all > > Just for fun, I have just downloaded the paper mentioned below and checked > it with R-1.6.1. > Everything is ok with exception of Table 2b, where I get always 1 instead of > 0.5: > >>pbinom(1e15,2e15,0.5) > > [1] 1 > > Which value should be correct? > > Best regards > Christian Stratowa > > =============================================> Christian Stratowa, PhD > Boehringer Ingelheim Austria > Dept NCE Lead Discovery - Bioinformatics > Dr. Boehringergasse 5-11 > A-1121 Vienna, Austria > Tel.: ++43-1-80105-2470 > Fax: ++43-1-80105-2683 > email: christian.stratowa at vie.boehringer-ingelheim.com > > > >>-------- Original Message -------- >>Subject: Re: [S] Exact p-values >>Date: Thu, 13 Feb 2003 18:31:38 +0100 >>From: "Rau, Roland" <Rau at demogr.mpg.de> >>To: 'Spencer Graves' <spencer.graves at PDF.COM>, Jose Mar?a Fedriani >>Laffitte <fedriani at ebd.csic.es> >>CC: s-news at lists.biostat.wustl.edu >> >>Dear all, >> >>in relation to your question, the following working paper of Leo Knuesel, >>University of Munich, might be of interest: >>"On the Accuracy of Statistical Distributions in S-Plus for Windows >>(1999)" >>You can download the paper from (pdf-Format, 45k): >>http://www.stat.uni-muenchen.de/~knuesel/elv/accuracy.html >> >>Best, >>Roland >> >> > -----Original Message----- >> > From: Spencer Graves [SMTP:spencer.graves at PDF.COM] >> > Sent: Thursday, February 13, 2003 6:12 PM >> > To: Jose Mar?a Fedriani Laffitte >> > Cc: s-news at lists.biostat.wustl.edu >> > Subject: Re: [S] Exact p-values >> > >> > >> > Try ( 1-pchisq(29.8, df=1)): With S-Plus 6.1, I got 4.78992e-008. >> > >> > By the way, the distribtion functions in R have more >>arguments. >> > For example, pchisq(29.8, df=1, lower.tail=F) produces the same >> > answer, and pchisq(29.8, df=1, lower.tail=F, log=T) produces its >>natural >> > logarithm. Also, pchisq, dchisq, qchisq, and rchisq in R all have an >> > "ncp" noncentrality parameter argument; only pchisq has such in S-Plus >> > 6.1. Similarly, none of the Student's t functions in S-Plus have a >> > non-centralitity parameter; in R, pt has an argument ncp, and from >>this >> > one can easily program ncp for dt, qt and rt. Also, the distribution >> > functions in the current release of S-Plus are known to have problems. >> > For example, pt(-1, Inf) = 0.5 in S-Plus 6.1, but 0.159 in R; >>clearly, >> > S-Plus gives a wrong answer without warning. >> > >> > Best Wishes, >> > Spencer Graves >> > >> > Jose Mar?a Fedriani Laffitte wrote: >> > >> > >Dear all, >> > > >> > > I want to get the exact p-values, on 1 degree of freedom, for an >> > array >> > >of chi-square values. When my chi-square values are equal or lower >>than >> > >29.7, I get the exact associated p-values. Thus, for instance: >> > > >> > > >> > > >> > >>pchisq(29.7, df=1) >> > >> >> > >> >> > >[1] 0.9999999 >> > > >> > >However, when my chi-square values are greater or equal to 29.8 what I >> > get >> > >is: >> > > >> > > >> > > >> > >>pchisq(29.8, df=1) >> > >> >> > >> >> > >[1] 1 >> > > >> > > >> > > Could anyone tell me how to fix this trivial issue? Very >>grateful, >> > Jose >> > >M. Fedriani >> > > >> > >**************************************** >> > >Jose M? Fedriani Laffitte >> > >Estacion Biologica de Donana (CSIC) >> > >Avda. M? Luisa s/n >> > >41013-Sevilla >> > >Spain >> > >Tel. +34-954232340 >> > >Fax +34-954621125 >> > >http://ebd.csic.es >> > > >> > >-------------------------------------------------------------------- >> > >This message was distributed by s-news at lists.biostat.wustl.edu. To >> > >...(s-news-request clipped)... > > >> > > >> > > >> > >> > >> > -------------------------------------------------------------------- >> > This message was distributed by s-news at lists.biostat.wustl.edu. To >> > ...(s-news-request clipped)... > > >>-------------------------------------------------------------------- >>This message was distributed by s-news at lists.biostat.wustl.edu. To >>...(s-news-request clipped)... > > >> > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > http://www.stat.math.ethz.ch/mailman/listinfo/r-help