[R] rbinom
Thomas Lumley
tlumley at u.washington.edu
Fri Feb 2 17:56:14 CET 2001
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list