[R] rbinom

[Thomas Lumley]

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
> 
>     if ( p == 0 ) p.is.zero <- p.is.zero + 1
>     if ( p == 1 ) p.is.one <- p.is.one + 1
> 
>   }
>   return(p.is.zero = p.is.zero,
>          p.is.one = p.is.one)
> }

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

	-thomas

Thomas Lumley			Asst. Professor, Biostatistics
tlumley at u.washington.edu	University of Washington, Seattle
^^^^^^^^^^^^^^^^^^^^^^^^ 
 NOTE NEW EMAIL ADDRESS

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