[R] binom.test inside ctest

[Peter Dalgaard BSA]

Emmanuel Paradis <paradis at isem.univ-montp2.fr> writes:

> At 00:22 12/09/00 -0400, you wrote:
> >I got 
> >
> >binom.test(65, 100, p=0.5)
> >
> >
> >p-value = 0.002654
> >
> >I thought in this case (when p=0.5 so a symmetric distribution)
> >we should get 
> >
> >p-value = 2*( 1-pbinom(64, 100, p=0.5) ) = 0.003517642
> >
> >Can somebody explain why the difference?
> >
> >Thanks.
> >
> >Mai Zhou
> 
> Hi,
...
> The relevant bit of code is:
> 
>                            ## Do
>                            ##   d <- dbinom(0 : n, n, p)
>                            ##   sum(d[d <= dbinom(x, n, p)])
>                            ## a bit more efficiently ...
>                            d <- dbinom(x, n, p)
>                            if(x / n < p) {
>                                i <- seq(from = x + 1, to = n)
>                                y <- sum(dbinom(i, n, p) <= d)
>                                pbinom(x, n, p) +
>                                    (1 - pbinom(n - y, n, p))
>                            } else {
>                                i <- seq(from = 0, to = x - 1)
>                                y <- sum(dbinom(i, n, p) <= d)
>                                pbinom(y - 1, n, p) +
>                                    (1 - pbinom(x - 1, n, p))
>                            }

Hum. Looks like we need some fuzz guards in there. The basic problem
is that

> sum(d[d <= dbinom(65, 100, .5)])
[1] 0.002653786
> sum(d[d <= dbinom(65, 100, .5)*(1+1e-7)])
[1] 0.003517642
> dbinom(35,100,.5)-dbinom(65,100,.5)
[1] 6.071532e-18

i.e. the "symmetrical" probability is ever so slightly bigger because
of inaccuracies, and so doesn't get counted in.

-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._