[R] rbinom

gb gb at stat.umu.se
Fri Feb 2 21:03:48 CET 2001


On Fri, 2 Feb 2001, Thomas Lumley wrote:

> On Fri, 2 Feb 2001, gb wrote:
>
> > Can someone tell me what I am doing wrong, or if there is
> > a bug in rbinom? I was simulating power in Fisher's
> > exact test and the normal approximation when I got strange
> > results. I have a stripped version here:
> >
> > -------------------------------------------------------
> > simulate <- function(replicates = 1000, size = 25,
> >                      p1 = 0.5, p2 = 0.5){
> >   p.is.zero <- 0
> >   p.is.one <- 0
> >   for (i in 1:replicates){
> >     x1 <- rbinom(1, size, p1)
> >     x2 <- rbinom(1, size, p2)
> >     p1 <- x1 / size
> >     p2 <- x2 / size
> >     p  <- ( p1 + p2 ) / 2
[snip]

> You should be getting p.is.zero about 1/4 of the time (and similarly for
> p.is.one), which agrees with what I get.
>
> Try plotting p1 & p2 against i.
>
> p1 & p2 are following a markov process  with absorbing states at 0 and 1,
> so they will quite quickly end up there.  If both end up at zero, you get
> p=0, if both at 1 you get p=1.
>
> My guess is that you didn't mean to call two different things 'p1'...

Thanks! Oh dear,... My only excuse is that it is late Friday night.

Thanks also to Matthew Austin for pointing out my error.

Göran

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._



More information about the R-help mailing list