[R] rbinom

[gb]

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)
}
----------------------------------------------------------
I run the function a few times and get
> simulate()
$p.is.zero
[1] 0

$p.is.one
[1] 0

> simulate()
$p.is.zero
[1] 0

$p.is.one
[1] 0

> simulate()
$p.is.zero
[1] 0

$p.is.one
[1] 0

> simulate()
$p.is.zero
[1] 975              ##!!!

$p.is.one
[1] 0

That is, I suddenly get highly extreme outcomes. And this behaviour
repeats. Sometimes with p == 1 a lot of times.

I run this on RH7.0, gcc 2.95.2, R-1.2.1

Göran
-- 
 Göran Broström                      tel: +46 90 786 5223
 professor                           fax: +46 90 786 6614
 Department of Statistics            http://www.stat.umu.se/egna/gb/
 Umeå University
 SE-90187 Umeå, Sweden             e-mail: gb at stat.umu.se

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