10 # pick from Urn 1 and put in Urn 2; is it blue?
11 if (runif(1) < nb1/n1) nb2 <- nb2 + 1
12 # pick from Urn 2; is it blue?
13 if (runif(1) < nb2/n2) count <- count + 1
14 }
15 return(count/nreps) # est. P(pick blue from Urn 2)
16 }Here is how we can do it without loops, usingapply():1 sim2 <- function(nreps) {
2 nb1 <- 10
3 nb2 <- 6
4 n1 <- 18
5 n2 <- 13
6 # pre-generate all our random numbers, one row per repetition
7 u <- matrix(c(runif(2*nreps)),nrow=nreps,ncol=2)
8 # define simfun for use in apply(); simulates one repetition
9 simfun <- function(rw) {
10 # rw ("row") is a pair of random numbers
11 # choose from Urn 1
12 if (rw[1] < nb1/n1) nb2 <- nb2 + 1
13 # choose from Urn 2, and return boolean on choosing blue
14 return (rw[2] < nb2/n2)
15 }
16 z <- apply(u,1,simfun)
17 # z is a vector of booleans but they can be treated as 1s, 0s
18 return(mean(z))
19 }Here, we set up a matrixuwith two columns ofU(0,1)random variates.
The first column is used for our simulation of drawing from urn 1, and the
second for drawing from urn 2. This way, we generate all our random num-
bers at once, which might save a bit of time, but the main point is to set up
for usingapply(). Toward that goal, our functionsimfun()works on one rep-
etition of the experiment—that is, one row ofu. We set up the call toapply()
to go through all of thenrepsrepetitions.
Note that since the functionsimfun()is declared withinsim2(), the locals
ofsim2()—n1,n2,nb1, andnb2—are available as globals ofsimfun(). Also, since
a Boolean vector will automatically be changed by R to 1s and 0s, we can find
the fraction ofTRUEvalues in the vector by simply callingmean().
Now, let’s compare performance.> system.time(print(sim1(100000)))
[1] 0.5086
user system elapsed
2.465 0.028 2.586
> system.time(print(sim2(10000)))310 Chapter 14