[R] FW: [Fwd: Re: [S] Exact p-values]

Fri Feb 14 13:35:04 CET 2003

Dear Spencer

Thank you for this extensive explanation of the problem.
I was just curious.

Best regards
Christian

==============================================
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-----
> From:	Spencer Graves [SMTP:spencer.graves at PDF.COM]
> Sent:	Friday, February 14, 2003 1:29 PM
> To:	Stratowa,Dr,Christian   FEX BIG-AT-V
> Cc:	r-help at stat.math.ethz.ch; David Smith
> Subject:	Re: [R] FW: [Fwd: Re: [S] Exact p-values]
> 
> 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
>